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1*.  AlSTKACT  (CsnMw*  CH  rtytfM  If  ntnmiy  tad  UmUtf  bf  Woe*  numbtf) 

The  prohlea  of  direction  finding  with  e  aenaor  errey  in  a  nultlple  source  envlroneent  Is 
addressed. 

In  this  resaarch  a  ganeralised  foraulation  to  the  problsn  is  developed.  In  the  generalized 
fomulation  a  aatrlz  pencil  is  eonstmcted.  The  rank  reducing  values  are  shown  to  contain 
Che  infomation  needed  to  estlnate  the  angles  of  arrival.  It  is  also  shown  chat  several 
different  nonsearch  procedures  are  special  cases  of  the  generalised  approach.  The  attractive 
ness  of  these  techniques  is  their  coaputatlonal  efficiency.  The  techniques  differ  depending 
upon  the  nanner  in  which  either  the  ■easursd  data  is  processed,  or  Che  angular  infomation 
is  eatracted  fro*  the  resulting  equations.  The  pencil  theoren  establishes  Che  relatlonahlp 
between  the  rank  reducing  values  and  the  functional  fom  generated  by  the  operators  applied 
to  Che  neasursnents.  This  generalised  approach  allows  for  the  Inclusion  of  prefiltering 
before  applying  an  operator  to  the  neasuranants. _ 
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CHAPTER  ONE 


li.TRODUCTION 

1.1  Problem  Stateaent 

The  problem  of  eaclmaclng  Che  angular  location  of  d  sources 
using  an  array  of  m  sensors  Is  addressed.  The  sources  are  assumed 
to  be  In  Che  far  field  with  respect  to  the  physical  size  of  the 
array.  The  signals  received  by  the  sensors  are  generated  by  the 
sources  themselves.  For  this  reason,  the  direction  finding  system 
used  In  this  work  Is  classified  as  passive.  The  general  system 
configuration  la  shown  in  Fig.  1.1. 

It  Is  assumed  that  the  model  which  governs  the  signal  received 
at  Che  sensor  Is  of  the  form 
d 

yj^(c,^)  •  I  ;  1-1,2,. ...m  (l.l-l) 

k-1 

where  Is  the  complex  envelope  of  the  signal,  £  Is  the  vector 

Q  —  02*****®jjJ 

where  the  1^  entry  6^  corresponds  to  the  angular  position  of  the 
1—  source,  Che  response  of  the  1—  sensor  to  the  1c^ 

source,  and  n^(c)  Is  Che  additive  noise.  This  noise  Is  the  sum  of 
external  and  Internal  noise.  External  noise  Is  assumed  to  be  received 
uniformly  from  all  directions  while  the  Internally  generated  thermal 
noise  Is  assumed  to  be  Identically  distributed  In  each  of  the  m 
channels. 
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Fig.  1.1  Passive  Direction  Finding  System 
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In  this  research  the  signals  generated  by  the  sources  are 
assuMd  to  be  narrowband.  A  signal  Is  classified  as  narrowband 
If  Its  bandwidth  Is  small  compared  to  the  Inverse  of  the  transit 
time  of  a  wavefront  across  the  array.  Given  N  snapshots  of  noise 
corrupted  data  collected  at  the  m  sensors >  the  problem  Is  to  estimate 
the  angular  position  of  the  d  sources. 

1.2  Previous  Results 

The  problem  of  estimating  the  bearing  or  direction  of  arrival 
(DOA)  of  radiating  sources  has  occupied  many  researchers  in  the  last 
two  decades.  This  problem  has  applications  In  areas  such  as  sonar, 
radar,  radio-astronomy,  and  seismology.  Also,  the  performance  of 
angle  of  arrival  estimators  plays  an  Important  role  In  determining 
the  ultimate  capability  of  any  surveillance  system.  Originally,  this 
problem  was  formulated  In  terms  of  the  classical  Bayesian  and/or 
Neyman-Fearson  decision  theories  [1],  [2].  However,  because  this 
approach  suffered  from  several  deficiencies,  a  variety  of  new  techniques 
have  been  proposed  In  recent  years.  In  this  section  we  review  some 
of  these  methods. 

Maximum  Likelihood  and  Least-Squares  Estimates 

Two  very  popular  approaches  for  constructing  parameter  estimates 
are  the  maxlnim  likelihood  and  the  least  squares  techniques .  Applica¬ 
tion  of  the  maximum  likelihood  approach  to  the  direction  finding  problem 
Is  discussed  by  Bohme  [3].  Let 

I*  [y^U,^),  y2(t,^),  ...,  y|^(t,£)l^ 
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T 

d«Tioc«  cooplcx  array  data  vector  where(*)  stands  for  the  transpose 
operator.  A  «■*.«■  likelihood  estiaate  (MLE)  can  be  found  provided  both 
the  probability  density  f(Y|9)  of  the  data  vector  7  is  known  and  the 
likelihood  function 


1(0)  •  log  f(Y[0) 

can  be  maximized  over  the  parameters 

[9^j  * 

A 

An  estiaate  ^  naxioizing  L(^)  over  £  is  defined  to  be  an  MLE  of 
Bbhoe  points  out  that  the  MLE  can  be  heavily  biased.  However,  for 
Independent  and  identically  distributed  random  variables  yj^Ct,^). 
y-(t,d) ,. . . ,y  (t,d)i  consistent  MLE  0  are  asymptotically  normally 
distributed  and  efficient.  This  means  that  the  distribution  of 
(£  •  6)  approaches  the  normal  distribution  with  zero  oean  and 
covariance  natrlx  given  by  the  inverse  of  the  Fisher  information 
matrix 

a  4nf(y, |6)  3  lnf(y.|£)  H 

J  -  '9  « - 55-^  )  <^-5? - > 

n 

where  (•}  stands  for  the  Hemltlan  operator.  Eg  means  expectation  with 
respect  to  6  and  3g/3£  1*  column  vector  with  elements  . 

An  estimate  of  ^  can  also  be  obtained  using  the  least,  squares 
approach.  For  the  narrowband  case,  the  data  vector  Y  has  the 
following  form 

Y  -  A5  +  N  (1.2-1) 


k 


w 


where  the  (m  x  d)  matrix  A  la  the  direction  matrix  whose  columns 
{a(6|^)t  k  ■  l»2«...,d}  are  the  signal  direction  vectors  for  the  d 
wavefronts . 

T 

S^[Sj^(c)f  s,(c}f...t  s^(t)]  where  Is  the  complex  envelope  of 

th 

the  le —  signal.  The  (m  x  i)  vector  M  represents  the  additive  noise 
at  the  m  sensors.  In  this  approach  one  attempts  to  fit  a  signal 
model  to  the  data  vector  ^  in  the  least  squares  sense;  that  Is, 
the  error 

E4  I  1  I  -  A  ll  ^ 

is  minimized  by  choosing  ^  and  ^  [4] .  If  we  set  the  derivatives 
of  £  with  respect  to  the  unknown  parameters  to  zero  in  trying  to 
minimize  E.  we  realize  that  some  of  the  necessary  conditions  are  in 
the  form  of  nonlinear  equations.  Attempting  to  solve  this  nonlinear 
set  of  equations  is  computationally  very  involved.  A  different 
approach  to  this  minimization  problem  was  suggested  by  several 
authors  [S  -  8].  The  idea  is  to  choose  values  of  ^  based  on  some  a 
priori  knowledge  or  based  on  some  preprocessing  and  find  the  associated 
^  that  minimize  E.  Then  the  values  of  £  can  be  altered  to  find  new 
values  for  ^  that  give  a  lower  mininum  for  E.  This  process  is  repeated 
until  a  local  mininum  for  E  is  found,  and  the  corresponding  ^  and  ^ 
are  the  estimates  of  the  unknown  parameters.  This  is  essentially 
a  search  procedure.  H.  Vang  [9]  comments  that  if  no  a  priori 
information  is  available,  the  local  minimum  that  is  first  reached  may 
give  an  estimate  which  is  far  from  the  true  values.  Both  approaches 
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discussed  above  prove  to  be  cooputatlonally  very  iavolved  casks. 


Beanforming 

Beanforming  Is  one  of  the  oldest  Ideas  In  array  processing 

for  determining  the  bearing  of  a  target.  A  simplified  block  diagram 

of  a  delay-and-sum  beamformer  la  shown  In  Figure  1.2.  The  Idea  Is  to 

align  Che  propagation  delays  of  a  signal  presumed  to  be  propagating 

In  a  direction  -k  so  as  to  reinforce  the  signal  [10].  For  example.  If 

the  sensor  delay  T^;  1*1, 2,..., m.  Is  Ideally  adjusted  to  compensate 

for  Che  signal  delay  *  k)/c,  where  ^  Is  the  spatial  location  of 

the  1—  sensor  and  c  is  the  speed  of  propagation,  the  signal  power 
2 

In  the  beam  Is  m  times  that  measured  at  each  sensor  while  the  beam 
iU3lae  power  Is  Increased  by  only  a  factor  of  m  (assuming  the  noise  to 
be  uncorrelated  from  sensor  to  sensor).  Dudgeon  [11]  explains  Chat 
Che  reason  for  studying  the  formation  of  beams  from  an  array  of  sensors 
Is  CO  use  the  signals  received  by  the  sensors  In  a  phased  manner  so 
as  to  preferentially  detect  signals  coming  from  a  particular  direction. 
In  addition,  by  averaging  over  many  sensors,  the  signal  to  noise  ratio 
(SNR)  Is  Increased  to  aid  In  Che  measurement  of  signal  parameters  [11]. 

The  delay-and-sum  beamformer  consists  of  computing  the  energy 
In  the  beam  for  many  directions  of  look  by  manipulating  the  delays. 
Maxima  of  this  energy  as  a  function  of  k  will  correspond  to  the  location 
of  Che  sources.  Observe  that  this  procedure  requires  computation  of 
the  energy  In  the  beam  at  every  direction  of  look.  For  this  reason, 
this  approach  can  be  classified  as  a  search  procedure.  The  energy 


la  th«  bean  when  steered  in  dlrecclon  k  is  given  by 

P(k)  -  A®  R  A  (1.2-2) 

D 

where  R  4E[YY  ]  is  the  spatial  correlation  matrix  of  the  sensor  outputs 
and  A  Is  the  direction  of  look  vector  with  elements  A^  ••  exp[J  ~  k*^] 
%>here  X  la  Che  signal  wavelength.  Assuming  the  noise  to  be  spatially 
white,  the  power  In  the  beam  when  steered  coward  the  source  is 

P(k)  -  m^o  ^  +  mo  ^  (1.2-3) 

“  3  n 

2  2 

where  and  o^  denotes  the  power  levels  of  the  signal  and  noise. 

The  value  of  k  which  yields  the  beam  power  given  In  (1.2-3)  Is  also 
known  as  the  Bartlett  esclnate  of  k.  For  this  value  of  k,  P(k)  Is 
maximized. 

The  resolution  of  this  approach  la  determined  essentially  by  the 
beam  pattern  of  the  array  of  sensors.  A  typical  beam  paccem  for 
a  delay  -  and  -  sum  beamfonaer  Is  shown  In  Fig.  1.3.  Using  this  classi¬ 
cal  method.  Increased  bearing  estimation  accuracy  can  only  be  obtained 
by  Increasing  the  aperture  of  the  array.  This  solution  la  of  limited 
utility  as  it  means  Increasing  the  physical  size  of  the  array. 

Don  H.  Johnson  (10]  observed  that  a  source  having  a  well  defined  bear- 
appears  to  be  coming  from  a  dominant  but  diffused  direction  as  well 
as  from  false  directions  covresponding  to  sldclobes.  The  sldelobes 
are  due  to  equal  weighting  assumed  for  each  aensor  output.  This 
approach  was  shown  to  be  Incapable  of  resolving  bearings  spaced  more 
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cloMly  than  th«  Baylclgh  Halt  regardless  of  SNR  [12].  This  alnlaum 
source  separation  below  which  sources  cannot  be  resolved  Increases  as 
the  magnitude  of  the  coherence  between  the  sources  Increases. 

Adaptive  Beamforalng 

A  way  to  alter  the  performance  of  the  benaformer  Is  to  weight  the 
sensor  signals  Individually  before  summing  them  as  shown  In  Figure  1.4. 

To  achieve  higher  resolution  modem  spectral  analysis  algorithms  are 
used  which  determine  the  weights  by  utilizing  the  data  measured  at 
the  sensors.  For  this  reason  such  systems  are  referred  to  as  adaptive 
beamformers.  The  problem  of  detenalnlng  the  sensor  weights  so  chat  the 
beam  pattern  has  some  desired  characteristic  Is  the  same  as  designing 
a  good  data  window  for  spectral  estimation »  or  designing  a  prototype 
low-pass  filter  for  use  in  a  digit*!  filter  bank  [11]  •  The  sensor 
weights  can  be  adjusted  to  maximize  the  SNR.  This  Is  analogous  to 
designing  a  Wiener  filter  given  the  spectral  estimates  of  the  signal 
spectrum  and  noise  spectrum  [II].  J.  P.Burg  [13]  notices  that  this 
filter  will  have  a  frequency  response  which  passes  chose  parts  of  the 
spectrum  where  the  SNR  Is  high  and  rejects  those  parts  where  the  SNR 
Is  poor.  This  adaptive  approach  performs  essentially  the  saise  as  the 
delay-and-sum  beamformer  when  the  source  angle  separation  Is  less  chan 
Che  besawldch.  Although  this  approach  has  better  resolution  capabilities 
chan  the  conventional  Fourier  transform.  It  Is  still  considered  poor. 

Spectral  Estimation 

The  propagating  waves  are  assumed  to  generate  »  homogeneous  random 
field  (l.e.,  a  field  that  is  space  stationary).  In  this  case  a 
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spectral  rapresantatlon  £or  the  field  exists •  similar  to  that  for 
stationary  random  processes,  which  consists  of  a  superposition  of 
traveling  waves  [14].  This,  In  effect,  explains  the  equivalence 
between  the  problem  of  determining  the  bearing  of  a  radiating 
source  with  an  array  of  sensors  and  the  problem  of  estimating  the 
spectrum  of  a  signal.  D.  Johnson  gives  an  excellent  review  of  this 
equivalence  [10]. 

Linear  prediction  techniques  can  also  be  applied  to  the  direction 
of  arrival  (OOA)  estimation.  A  tutorial  review  of  linear  prediction 
is  given  by  J.  Makhoul  [15].  Ulrych  and  Clayton  [16]  as  well  as  Tufts 
and  Kumaresan  [17]  have  applied  the  forward-backward  linear  prediction 
method  to  the  narrowband  linear  array  case.  Among  the  most  well  known 
high  resolution  spectral  estimation  methods  are  Capon's  maximum  likeli¬ 
hood,  Burg's  maximum  entropy  and  the  linear  predictive  method.  Ve 
shall  give  a  brief  discussion  of  these  three  techniques  as  applied  to 
DOA  estimation. 

Capon's  maximum  likelihood  method  does  not  make  use  of  the 
standard  MLE.  Instead,  a  constrained  optimization  problem  Is  solved. 

The  direction  of  look  A  is  found  idilch  yields  the  beam 

energy  P(fc)  given  In  (1.2-2)  subject  the  the  constraint  A^  ■  1,  where 
E  represents  an  Ideal  plane  wave  corresponding  to  the  direction  of  look. 
This  constraint  fixes  the  processing  gain  to  1  for  each  direction  of 
look.  Minimizing  the  beam  energy  tends  to  reduce  the  contributions 
to  this  energy  from  sources  and  noise  not  propagating  in  the  direction  of 
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look  [10].  For  Chls  reason  S.  R.  Degraaf  and  D.  H.  Johnson 
called  chls  approach  sdnlaun  energy  adaptive  beamfomlng .  Using 
Che  Lagranglan  approach*  Che  solution  of  Che  above  staced  optimiza¬ 
tion  problem  is  given  by 


A 


(1.2-4) 


where  R  is  Che  special  correlacion  matrix  of  the  data  vector  T 
which  was  defined  previously.  This  approach  is  classified  as  a 
search  procedure  since  the  vector  A  has  to  be  computed  for  each 
direction  of  look.  When  the  noise  is  assumed  to  be  white  and  Che 
beam  is  steered  toward  the  source,  the  beam  power  estimate  is 


P(k)  - 


2  2 

where  a  and  a  denote  the  power  levels  of  the  signal  and  noise. 

9 

Burg  [18]  suggested  using  the  principal  of  maximum  entropy  to 
estimate  the  DOA'a.  For  a  given  direction  of  look,  the  entropy 

H  -  j  in  FOc)  dk  (1.2-5) 

of  the  power  spectrum  P(k)  is  obtained.  Burg  has  shown  chat 
maximizing  the  entropy  is  equivalent  to  choosing  the  "most  likely" 
spectrum.  For  Che  case  of  an  equally  spaced  linear  array,  the 
naxlaum  entropy  solution  was  shown  to  be  equivalent  to  the  linear 
prediction  solution. 
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S.  V,  Lang  and  J.  H.  Hcllalan  [19].  among  others,  applied 
linear  prediction  (LP)  to  array  processing  problems.  The  idea  la 
to  estimate  the  output  of  the  sensor  es  a  linear  combination  of  the 
other  sensor  outputs 

•  -  I  (c.e).  (1.2-6) 

tti 

The  linear  predictive  coefficients  are  obtained  by  again  solving  a 
constrained  optimization  problem,  the  constraint  being  ~  1.  Once 
these  coefficients  are  known,  the  power  spectrum  P(k)  is  obtained. 

The  power  spectrum  P(k)  is  computed  for  different  directions  of  look 
and  Che  source  bearings  correspond  to  the  maxlmas.  Observe  chat 
this  la  also  a  search  procedure.  An  open  question  in  the  LP  approach 
Is  Che  choice  of  1,  l.e.  the  choice  of  the  sample  Co  be  predicted  by 
the  other  samples. 

The  three  methods  described  above  are  capable  of  resolving  closely 
spaced  source  bearings  if  the  SNR  is  large  enough.  Using  incoherent 
sources  and  a  linear  array  containing  10  sensors,  it  was  shown  that 
the  linear  prediction  processing  algorithm  is  uniformly  the  moat 
capable  of  resolving  closely-spaced  signals  [12].  Thus,  these  methods 
have  better  resolution  capabilities  than  the  beamformers  presented 
earlier.  However,  one  mejor  drawbeck  of  these  high  resolution  pro¬ 
cedures  is  they  are  computationally  very  complex  due  to  their 
search  nature.  Finally  it  should  be  noted  that  spectral  estimation 
methods  are  not  applicable  to  the  DOA  estimation 


probl**  «h«n  th«  sourca  signals  ars  corralatad  because  the  wave- 
field  generated  by  correlated  signals  Is  not  homogeneous. 

Slgnal-Subspece  Processing 

In  recent  years,  there  has  been  a  great  deal  of  Interest  In 
spectral  estimation  procedures  based  on  an  eigenvalue-eigenvector  de¬ 
composition  of  the  spatial  correlation  matrix.  In  the  mathematical 
literature,  this  Is  known  as  principal  component  analysis.  R.  Kumaresan 
[20]  points  out  that  principal  component  analysis  Involves  linearly  trans¬ 
forming  the  data  along  the  first  few  principal  eigenvectors  of  the 
covariance  matrix  of  the  data  [21].  The  original  work  In  this  area 
Is  due  to  H.  Hotelling  [22].  C.  R.  Rao  [23]  showed  chat  several 
problems  In  mulclvaclate  analysis  have  solutions  based  on  calculating 
the  first  few  eigenvectors  of  a  correlation  matrix.  In  array  pro¬ 
cessing,  Interest  stems  from  the  fact  Chat  there  are  relations  between 
Che  eigenvectors  of  the  spatial  correlation  matrix  and  the  directions 
of  arrival.  Pisarenko  [24]  was  the  first  to  realize  the  Important  prop¬ 
erties  of  Che  eigenvectors  of  the  correlation  matrix.  This  Idea  was 
Chen  vigorously  developed  by  Cantonl  and  Godora  [25],  Owsley  [26], 
Liggett  [27],  Reddl  [28],  Schmidt  [29]  and  Bienvenu  [30]. 

All  these  slgnal-subepeee  approaches  assume  the  background  noise  to 
be  Independent  from  sensor  to  sensor  end  the  spatial  coherence  matrix 
CO  be  Che  Identity  matrlz.  It  should  be  noted  that  when  the  noise 
is  nomdilce  but  has  known  covariance,  the  problem  can  still  be  handled 
through  prewhitenlng.  A  Faulraj  and  T.  Kailach  [31]  proposed  a 
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solution  to  th«  unknown  nola*  covnrlanc*  problea  when  the  noise 
field  Is  Invariant  under  different  array  positions  (l.e.,  the  noise 
covariance  matrix  Is  unchanged ) .  Situations  where  this  assumption  Is 
valid  are  not  uncononn  In  sonar  applications. 

Ue  now  present  the  building  block  of  a  standard  elgenstruccure 
algorithm.  Using  (1.2-1),  the  spatial  covariance  matrix  can  be 
written  as 

R  -  ECl  l“l  -  ASA®  +  a^l  (1.2-7) 

where  S  Is  the  covariance  matrix  of  the  signals  s(t),  l.e. 

S  -  E(s(t)  s®(t)I. 

2 

a  Is  an  unknown  constant  representing  the  noise  power  at  each  sensor, 

and  E(*)  denotes  expectation.  Let  (X,  >  X.  >. . .  >  X  }  and 

I  z  n 

e^,  ...,  be  the  eigenvalues  and  the  corresponding  eigenvectors 
of  R.  One  can  show  that  If  S  la  nonsingular  (l.e.  the  source  signals  are 
not  coherent) ,  and  If  the  number  of  sources  d  Is  less  chan  the  number  of 
sensors  m,  then 

1)  The  minimum  eigenvalue  of  R  Is  a  with  multiplicity  (m-d) ,  l.e. 

Vl  "  ^d+2  "  •  ••  “ 

2)  The  eigenvectors  correspondli^  Co  the  minimal  eigenvalues 
are  orthogonal  to  the  columns  of  Che  matrix  A  or  the  signal 
direction  vectors.  In  particular. 
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{ t  S^^2*****  ^^^  *^  *  ■  *  *  **^^  (1*2‘“9) 

where  the  symbol  j.  denotes  orthogonelity . 

The  key  steps  of  this  elgorlthm  ere  as  follows.  First  determine 
the  number  of  sources  d  using  (1.2-8).  Next  find  the  noise  subspace 
as  the  space  of  the  eigenvectors  corresponding  to  the  smallest 
eigenvalue  of  R.  Finally,  the  desired  direction  vectors  (hence  the 
OOA's)  are  found  by  searching  the  array  manifold  to  locate  those 
vectors  that  are  orthogonal  to  this  noise  subspace  [29].  Like  earlier 
approaches,  these  signal -subspace  techniques  are  also  search  procedures. 

The  above  algorithm  assumes  that  ve  have  perfectly  known  eigen¬ 
vectors  and  eigenvalues.  In  practice,  however,  R  has  to  be  estimated 
from  the  array  data  vector  Y.  Due  to  errors  In  estlx&atlng  R  from 
finite  samples  and  because  of  finite  precision  arithmetic,  the 
eigenvalues  and  eigenvectors  will  be  perturbed  from  their  true  values. 
There  Is,  therefore,  zero  probability  that  the  smallest  eigenvalue 
of  R  will  have  multiplicity  (m-d),  l.e.  the  small  eigenvalues  will  all 
be  different  with  probability  one.  A  more  sophisticated  approach,  based 
on  statistical  considerations.  Is  needed  to  determine  the  underlying 
multiplicity  [32].  The  problem  of  estimating  the  number  of  sources 
present  Is  an  on  going  one.  Hany  solutions  have  been  proposed  but 
none  Is  fool  proof.  Liggett  [27]  was  the  first  researcher  to  use 
eigenvalues  of  the  estimated  R  to  determine  the  number  of  signal  sources. 
Be  fitted  a  multivariate  signal  model  to  the  observed  data.  Bartlett 
[33]  and  Lawley  [34]  developed  a  procedure  based  on  a  nested  sequence 
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of  hypothesis  tests.  For  each  test  the  likelihood  ratio  statistic 
Is  covputed  and  compared  to  a  threshold.  The  hypothesis  accepted 
Is  the  first  one  for  which  the  threshold  Is  crossed.  The  problem 
with  this  method  Is  the  subjective  Judgment  required  for  deciding  on 
the  threshold  levels  (32].  One  of  the  most  celebrated  approaches 
to  this  problem  Is  the  Alkake  Information  theoretic  criteria  (AlC) . 
The  advantage  of  this  approach  la  that  no  objective  judgment  Is 
required  In  the  decision  process  [35].  The  minimum  description 
length  (MDL)  Is  another  approach  that  does  not  require  objective 
judgment.  Recently,  a  new  approach  based  on  the  AIC  and  MDL  has  been 
proposed  by  M.  Wax  and  T.  Kallath  [32].  They  view  the  problem  as  a 
model  selection  problem  and  then  apply  the  AIC  or  MDL  for  model 
selection. 

To  close  this  subsection  we  present  one  of  the  most  promising 
subapace  approaches ,  called  Multiple  Signal  Classification  (MUSIC) , 
which  was  proposed  by  Schmidt  [36].  MUSIC  can  be  shown  to  provide 
asymptotically  unbiased  estimates  of 

1)  number  of  Incident  wavefronts  present 

2)  directions  of  arrival 

3)  strengths  and  cross  correlations  among  the  incident  waveforms. 

4)  nolse/interference  strength 

5)  polarizations 

Given  the  m  x  a  spatial  covariance  matrix  R  of  the  data  vector  Y, 
Schmidt  showed  that  the  eigenvectors  associated  with  the  minimum 
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eigenvalues  of  R  are  orthogonal  to  the  space  spanned  by  the  coluons 
of  A,  l.e.  the  signal  direction  vectors.  Thus,  Schmidt  proposed 
searching  In  the  a^(6)  continuum  for  the  values  of  a(9)  which  will 
maximize  the  expression 

P(e)  - -5 - - -  (1.2-10) 

a“(9)E^Ej  a(e) 

where  Is  defined  to  be  the  a  x  (m-d)  matrix  whose  columns  are  the 
(a-d)  noise  eigenvectors.  These  new  eigenvectors  are  the  eigenvectors 
corresponding  to  the  (m-d)  smallest  eigenvalues.  Once  the  signal  direction 
vectors  are  known,  the  OOA's  can  then  be  obtained.  Note  that  MUSIC 
necessitates  searching.  Schmidt  generalized  this  approach  to  Include 
the  added  complexity  of  signal  polarization. 

Schmidt  compared  MUSIC  with  the  conventional  beaztformer,  Capon's 
maxloum  likelihood  and  Burg's  maximum  entropy  approaches  (See  Fig.  1.5). 
MUSIC  was  found  to  have  superior  bias,  error  variance  and  resolution 
performance.  However,  like  all  the  earlier  search  methods,  this  super- 
resolution  approach  la  very  Inefficient  computationally. 

Nonaearch  Procedures 

Currently,  nonaearch  procedures  are  being  developed  [37],  [38], 

[39].  These  algorlthaw  seem  to  have  important  advantages  over  search 
procedures.  Aa»ng  these  advantages  are  [37]: 

1)  They  are  computationally  leas  complex  because  a  search  procedure 
Is  not  needed. 

2)  These  algorithms  do  not  require  knowledge  of  element  characteristics 
(l.e.  directional  pattern,  gain/phase  response). 
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Fig.  1.5  ExAapld  of  Azimuth  Only  Dir«ctlon  Finding  Performance. 
(From  [291) . 
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3)  They  do  not  require  a  calibration  of  the  array.  This 

conpletely  ellalnaces  need  for  storage  of  the  array  manifold 
which  be  very  large  for  multidimensional  problems . 

This  dissertation  concerns  itself  with  these  newly  developed 
nonsearch  procedures.  Although  the  techniques  are  developed  for 
narrowband  signals t  wideband  signals  can  also  be  handled  by  decompos' 
ing  them  into  narrowband  signal  sets  using  comb  filters.  A  discussion 
of  OOA  estimation  of  wideband  sources  is  given  by  H.  Wang  and  M.  Kaveh 
[401. 


1.3  tesearch 

Monaearch  methods  for  solving  the  problem  of  direction  finding 
with  a  sensor  array  in  a  multiple  source  environment  are  the  main 
topic  of  this  research 

Chapter  2  is  devoted  to  formulating  a  generalized  approach  for 
these  nonaearch  procedures.  The  generalized  approach  consists  of 
applying  an  operator  to  the  received  signals  in  order  to  form  a 
matrix  pencil  M-IN.  The  rank  reducing  values  of  A.  are  shown  to  con¬ 
tain  the  information  needed  to  estimate  the  angles  of  arrival.  The 
pencil  theorem  proposed  in  this  chapter  establishes  the  relationship 
between  the  rank  reducing  values  of  X  and  the  functional  form  f((^^) 
generated  by  the  operators  applied  to  the  measurements.  When  the 
matrices  M  and  N  are  square,  the  rank  reducing  values  of  X  are  the 
generalized  eigenvalues  of  the  pencil  M-XN.  Two  methods  are  proposed 
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to  conpuce  these  values  of  X  when  M  and  N  are  nonsquare. 

Several  different  nonsearch  procedures  are  special  cases  of  the 
generalized  approach.  These  methods  differ  primarily  depending  upon 
the  operator  used.  In  chapter  3  two  different  techniques  are  pre¬ 
sented  and  fonulated  In  terms  of  the  generalized  framework.  The 
first (  ESPRIT •  makes  use  of  a  phase  delay  operator  while  the  second 
makes  use  of  a  summation  operator. 

In  chapter  4  the  moving  window  operator  Is  analyzed  In  detail. 

This  operator  Is  also  formulated  within  the  generalized  framework 
and  is  shown  to  hold  for  coherent  sources.  A^t  first  the  moving 
window  Is  presented  for  the  case  of  deterministic  signals.  In 
section  4.1-2  the  moving  window  Is  shown  to  be  applicable  to  the 
zero  mean  random  signal  case.  In  section  4.2  the  moving  window  oper¬ 
ator  Is  applied  to  a  rectangular  planar  array.  As  before »  the  rank 
reducing  numbers  of  the  matrix  pencil  N-XN  contain  the  Information 
needed  to  locate  the  sources.  Sections  4.3  and  4.4  present  Prony's 
and  Pisarenko’s  algorithms.  Their  relationship  to  the  moving  window 
Is  developed.  In  section  4.5  computer  simulation  results  are  pre¬ 
sented.  The  estimation  accuracy  of  the  Prony,  Pisarenko,  and  moving 
window  methods  are  evaluated  In  terms  of  bias  sad  variance  based  on 
several  Monte  Carlo  runs.  The  moving  window  is  also  compared  to  ESPRIT. 

Finally,  a  summary  and  suggestions  for  future  work  are  given 
In  chapter  5. 
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CHAPTER  2 


THE  GENERALIZED  APPROACH 

As  dlscusssd  In  ths  Introduction,  several  high  resolution 
spatial  processing  algorithms  such  as  the  likelihood  method, 

maximum  entropy  method,  linear  predictive  method,  and  MUSIC  have  been 
proposed  for  estimating  the  directions  of  arrival  of  Incoming  signals. 
These  techniques  are  essentially  search  procedures  and,  because  of 
chat,  are  computationally  very  complex.  Honsearch  procedures  have 
recently  been  proposed  to  alleviate  this  computational  burden. 

In  this  chapter  a  generalized  formuatlon  Is  proposed  for  these 
newly  developed  techniques.  The  pencil  theorem  which  forms  Che 
foundation  of  this  work  Is  developed  In  section  2.1.  This  theorem 
states  Chat  the  rank  reducing  values  of  a  certain  matrix  pencil  con¬ 
tain  the  Information  needed  to  estimate  the  angular  location  of  Che 
d  sources.  Section  2.1  discusses  the  procedures  one  can  use  to 
obtain  these  rank  reducing  numbers. 

2 . 1  The  Pencil  Theorem 

The  basic  problem  under  consideration  Is  that  of  estimating  the 
angular  location  of  d  sources  given  measurements  from  an  array  of  m 
sensors.  The  measurcuents  are  modeled  by  a  linear  combination  of  d 
exponentials  whose  exponents  J^^,l  ■  1,2,...,  d  contain  the  Informa¬ 
tion  about  Che  angular  location  of  the  sources.  The  generalized 
approach  Is  based  on  the  concept  of  a  matrix  pencil  which  Is  defined 
as  follows  [41].  Let  M  and  M  be  two  k  x  1  matrices.  Also,  denote 
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Che  epece  of  cooplex  miobere  by  The  sec  of  ell  mecrlces  of  the 
form  M  -  XN  irlch  Xe^  ie  seld  Co  be  e  pencil. 

The  mecrlces  M  end  N  ere  required  co  have  che  decoiqraslclons: 

M  -  EF  end  N  -  EDF  (2.1-1) 

where  E  Is  e  k  x  d  mecrlz  (k  ^  d),  F  Is  e  d  x  I  mecrix  (1  ^  d), 

D  Is  e  d  X  d  dlegonel  aacrlx,  end  E.  D,  end  F  ere  ell  of  renk  d. 

In  eddlclon,  che  11—  encry  of  che  dlegonel  mecrix  D  Is  required  co 
be  of  che  form 

dfi  -  -  f[g(e^)];  1  -  1,2, ...d  (2.1-2) 

where  6^  Is  che  enguler  loceclon  of  che  1—  source.  The  mecrlces 
M  and  N  end  che  funcclon  f(*)  ere  decermlned  by  che  operecor  T{*} 
applied  CO  Che  signals  received  ec  che  m  sensors.  In  general,  dlfferenc 
approaches  may  employ  dlfferenc  operecors  T{*}. 

The  following  cheorem  escebllshes  che  releclonshlp  becween  che 
values  of  X  for  which  che  rank  of  che  pencil  M  -  XM  decreases  by  1 
and  che  funcclons  f  (*)  which  conceln  Che  dlrecclons  of  arrival  (DOA) 
Informaclon. 

Pencil  Theorem; 

Lee  M  -  XN  be  a  macrlx  pencil  where  M  and  N  have  che  decomposlclons 

N  •  EF  and  N  •  EDF 

as  deed  above.  Then,  che  values  of  X  which  decrease  che  rank  of  che 
pencil  M  -  XN  by  1  are  given  by 
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■r- 


-  {f[g(e^)ir^  ;  1-1 . d. 


Proof; 

Glvan  ch«  aatrlx  ponell  M  -  XM.  It  follow*  chat 

M  -  XN  -  EF  -  XEDF  -  E(I  -  XD)F  (2.1-3) 

uh«r*  th«  11—  entry  of  the  dlagonel  matrix  (I-XD)  Is  given  by 

C(I  -  ^ 

If  X  ■  (f(<t>^)}~^«  then  the  1—  column  of  (I  -  XD)  becomes  zero.  The  rank 
of  (X  -  XD) t  denoted  by  rCl  -  XD],  Is  thus  reduced  by  1  to  (d-1) .  In 

general 

r(M-XNl  -  rCE(I  -  XD)Fl  -  min  {r[EJ,  r[I-  XD] ,  r[r]}. 

However,  by  assumption 


riE]  -  r[F]  -  d. 

Therefore,  the  rank  of  the  pencil  M  -  XN  Is  decreased  by  I  when 

\  -  {f(^^)r^  ;  1  -  1,2, ...d.  (2.1-4) 

2.2  Evaluation  of  the  Rank  Seducing  Numbers 

The  procedure  used  to  compute  the  X^'s  Is  now  discussed.  It  Is 
shown  to  depend  on  whether  M  end  N  ere  square  or  non-square  matrices. 
Denote  the  determinant  of  a  matrix  M  by  det(M),  and  the  space  of 
complex  numbers  by  f.  In  ease  M  end  N  arc  SQUARE  matrices,  the 
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generalized  eigenvalues  of  che  sacrlx  pencil  M:1  '  are  idencical 


CO  chose  obtained  In  (2.1-4)  for  Che  matrix  pencil  K-IS  when  fl  and  N 
are  square. 

2.  A.  Granoaian  approach  can  also  be  used  Co  determine  che  a^'s.  Denote 
che  i—  column  of  M  -  IM  by  (M  -  XN)^.  Also,  define  che  inner  product 
of  two  vectors  X  and  ^  by 

<X,  1  >  -  (2.1-3) 

This  approach  consists  of  checking  che  dependence/independence  of  che 
columns  of  che  matrix  pencil  M  -  XN.  To  this  end  a  Gran  matrix  G 
wnose  ij —  entry  is  given  by 

G^j  -  <(M  -  aS)^,  (M  -  .\S)^> 

is  formed.  Conpucing  che  deceminanc  of  G  results  in  a  polynomial  P(a) 
whose  non-zero  zeros  are  che  rank-reducing  numbers  of  M  -  XN. 

Having  shown  how  che  rank  reducing  numbers  or  a  matrix  pencil  can 
be  decemined,  we  now  show  chac  different  nonsearch  procedures  for  DOA 
esciaacion  can  be  foroulacad  within  this  generalized  framework. 
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CHAPTER  3 


VARIOUS  OPERATORS  FOR  USE  IN 
THE  GENERALIZED  APPROACH 

Th*  rank  reducing  numbers  of  Che  aacrlx  pencil  M  -  XN  will  be 
functions  of  the  angular  locations  of  the  d  sources  provided  the 
1~  diagonal  clement  of  the  matrix  0  has  the  functional  form  given 
by  (2.1-2).  This  functional  form*  however.  Is  not  dictated  by  the 
location  of  the  sources  but  by  the  operator  used  to  process  the  data 
received  at  the  sensors.  Different  methods  formulated  within  the 
framework  of  the  generalized  approach  use  different  opracors. 

In  this  chapter  two  direction  of  arrival  csclmaclon  techniques 
are  presented.  The  first,  ESPRIT,  which  makes  use  of  a  delay  oper¬ 
ator,  la  discussed  In  section  3.1.  The  second  approach,  presented 
In  section  3.2,  is  a  generalization  of  a  system  Identification  scheme 
originally  proposed  by  Jain  [42].  This  scheme  uses  a  summation  oper¬ 
ator. 


3.1  Phase  Delay  Operator  (ESPRIT) 

ESPRIT  [37]  la  a  subapace  approach  to  direction  of  arrival  estima¬ 
tion  which  employs  a  planar  array  of  arbitrary  geometry  composed  of  m 
matched  sensor  doublets  whose  elements  are  tranalatlonally  separated 
by  a  known  constant  displacement  vector  A  (See  Fig.  3.1).  In  terms 
of  the  generalized  formulation.  It  will  be  shown  that  ESPRIT  uses  a 
phase  delay  operator  which  results  In  the  function 

f(d^)  -  e  (3.1-1) 
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Xhtt  array  of  sanaors  la  convanlancly  daacribad  aa  being  coaq)rlaad  of 
two  aubarraya  X  and  Y,  Identical  in  every  aapect  although  phyalcally 
dlaplaced  from  each  other  by  a  known  dlaplacencnt  vector  A.  Asatiae 
there  are  d  narrowband  atatlonary  aourcea  located  at  aaimuthal  anglea 
6^,  1-1. 2,... d.  which  are  impinging  on  the  array  aa  planar  wavefronta 
and  emitting  aignala  whoae  compltft  envelopea  are  denoted  by  s^(t) ,  k- 

l.  2.....d.  To  be  able  to  aolve  for  the  DOA'a  d  haa  to  be  leaa  than 

m.  The  aignal  received  by  the  i^  aenaor  ia  the  superpoaltlon  of  d 


impinging  wavefronts  plus  zero-mean  additive  noise.  Define  e^(d^)  to 

be  the  relative  response  of  the  1~  sensor  at  either  subarray  to  the 

source,  to  be  the  center  frequency  of  each  of  the  spatial 

sources,  c  to  be  the  speed  of  propagation  of  the  plane  waves,  and  n  (t) 

ch  ^ 

and  n  (t)  to  be  the  additive  noise  at  the  elements  In  the  1—  doub- 
^1 

let.  Then  the  output  of  the  1~  sensor  can  be  expressed  as. 


X.  (t,9)  -  I  s,,(t)  +  o_  (t) 

^  k-l  ^  1 


i 

y  (t,9)  -  y  s,(t)  (t) 

1  k-l  'i 


where 

9  - 

(® j »  ®2* * • • 

and 

“o 

"  Q  ^  ®k' 

(3.1-2) 

Let 

x'^(t)  -  [x^(t,9),  x^U.i) . 

i^(t)  -  [yj(t.9),  y^U.e).....  y,(t.9)l 
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b«  sinulcaocously  saopled  daca  vaccors  of  both  subarrays.  Similarly » 
dafina 

S^(t)  -  [8^(t),  s^Ct),....  aj(t)] 

n^(t)  -  [n  (t),  n  (t),...,  n  (t)] 

*1  *2  *« 

a^(t)  -  [Uy  (t),  Uy  (t)]. 

Dafina  A  to  ba  tha  m  x  d  dlractlon  matrix  whose  columns  are  the  signal 
direction  vectors  for  the  d  wavefronts  defined  by 

a^  (6^)  ■  [Sj^  (6j^)  1*2  (0j^)  t  *  •  • »  ^ »  k“l»...d. 

Also,  define  to  ba  the  diagonal  matrix  of  the  phase  delays  between 
the  doublet  sensors  given  by 

J«,  J«2 

♦  ■  dlag  [e  *■,  e  e  “] .  (3.1-3) 

It  follows  that 

X(t)  -  A  S  (t)  +  n^(t) 

Y(t)  -  A  ♦  S(t)  +  n  (t).  (3.1-A) 

_  -  “y 

As  a  first  step  in  this  algorithm  the  minimum  description  length 
(MDL)  criterion  [32]  is  used  to  estimate  d,  the  number  of  sources. 

Using  both  this  estimate  and  tha  assumption  that  the  noise  is  white 
and  uncorrelated  from  sensor  to  sensor,  the  maximum  likelihood  estimate 
of  the  noise  power  is  obtained.  The  matrix 


31 


pencil  C-XB  is  then  formed  uhere  C  is  the  autocorreletion  of  the  data 
received  by  subarray  X  and  B  is  the  cros8<-corrclation  matrix  between 
the  data  received  at  subarrays  X  and  Y.  Note  that  C  and  B  are 


square  correlation  matrices.  Aasuning  that  £(t) .  n^(t)  and  n  (t)  are 

— y 

statistically  independent •  then 

C  -  EtX(t)  X®  (t) ) 

-  £((A  S(t)  +a  (t))(AS  (t)  +n  (t))®] 

-  A  S  A^  +  R  (0)  (3.1-5) 

n  Q 
-Tt-Jt 

where 

S  -  E[S(t)  S®  (t)]  (3.1-6) 


is  the  d  X  d  correlation  matrix  of  the  signals  s(t)> 


K  n  -  eco  (c)  (t)] 

n  n  —x  —X 


la  the  d  X  d  correlation  matrix  of  the  noise  and  E(’)  denotes  the 
expected  value. 

Let 

E  -  AS  and  F  -  A^. 


Using  these  expressions  C  can  be  rewritten  as 


C  -  EF  R  „  (0) 

n  n 

-X— X 


(3.1-7) 


Define 


0  •  ♦  ■  dlag  [e  »  e  .....e  “) .  (3.1-8) 
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A 


Observe  that  the  matrix  pencil  formed  in  this  manner  has  the  decom¬ 
position  required  by  the  generalized  formulation. 

Recall  that  the  m  x  d  matrix  E  is  the  product  of  the  direction 
matrix  A  and  the  signal  covariance  matrix  S.  The  matrix  E  Is  required 
to  be  of  full  rank.  This,  In  turn,  requires  that  the  matrices  A  and 
S  both  be  of  rank  d.  The  m  x  d  matrix  A  Is  of  rank  d  as  long  as  the 
signal  angles  of  arrival  are  distinct.  However,  the  rank  of  S  will 
be  less  than  d  when  the  signals  are  coherent.  This  explains  why 
ESPRIT  falls  in  the  case  of  coherent  signals.  The  matrix  F  Is  of  full 
rank,  as  required,  since  It  Is  the  matrix  A  .  The  diagonal  matrix 
0  given  In  (3.1-8)  is  also  of  rank  d. 

As  mentioned  In  (3.1-1)  the  phase  delay  operator  of  ESPRIT  results 
In 

f((ti^)  •  e 

From  the  pencil  theorem  It  follows  that  the  G.E.'s  of  M-XN  are  related 
to  by 

Xi  -  U(ii^)}  -  e  1  •  1,2,. ..,d. 

Once  Che  G.E.'s  are  known,  the  DOA's  can  be  obtained  from  (3.1-2)  using 
Che  relation 

9.  •  arcsin  {-J  InA  }. 

1  ^ 

Mote  Chat  the  argument  of  arcsin  is  real  since  in  is  purely  imaginary. 

It  can  be  concluded  chat  this  nonscarch  procedure  is  one  application 
of  Che  generalized  approach.  The  operator  used  by  ESPRIT  is  the  phase 


delay  operator  which  results  in  the  functional  relationship  given 
in  (3.1-1).  ESPRIT  has  all  the  advantages  discussed  in  the  intro¬ 
duction,  i.e.,  speed,  storage,  and  indifference  to  calibration. 


3.2  S**™“»tion  Operator 

In  this  section  ve  generalize  a  system  identification  scheme  first 
proposed  by  V.  R.  Jain  [42].  In  this  generalization  a  linear  array 
composed  of  m  identical  omnidirectional  sensors  with  uniform  spacing 
A  la  used  to  estimate  the  direction  of  arrival  of  far  field  point 
sources.  Assume  there  are  d  ^  m/2  narrowband  deterministic  sources 
with  spectra  centered  at  (Uq  and  located  far  enough  from  the  array 
such  that  their  wavefronts  impinging  cn  the  array  are  planar.  ®k 
the  direction  of  arrival  of  the  source  and  a^(c)  is  the  additive  noise 
at  the  i^  sensor.  Additive  noise  is  present  at  all  m  sensors.  The  signal 
received  at  the  1^^  sensor  is  expressed  as 

<i)- 

d  J—  A{l-l)3in6, 

y^^Ct,^)  ■  I  ^  +Qj^(t);  1-1, 2,..., m 

k-1 


As  before,  denotes  the  complex  envelope  of  the  signal  emitted 

1.  u  1  th 

by  the  source. 

For  a  reason  which  will  become  apparent  later  the  i—  measurement 
y^(t,£}  is  weighted  by  a  known  decaying  exponential  e  The 

weighted  signal  is  given  by 


u)^(t,9^) 


-bvl-l) 


yj^(t,9) 


**  i\  J  7?- A(i-l)  8ln0.  . 

-  I  s^(t)  e-'»<^-^>  e  C  '^  +  e-*’<^-'>n,(t). 


k-1 
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Let 

♦k  *  t-b  +  J  ainej^]  (3.2-1) 

Hence,  (i>^(t,^)  cen  be  expressed  as 

a^(t.8)  -  8j^(t)(g(e^)]^^“^^  +  ,-b(i-l) 

Aaaumtng  Che  signals  Co  be  decermlnlsclc  and  the  noise  Co  be  zero 
mean>x^(c,^)  Is  defined  as 

.^(t.9)  -  Et»^(t.e,l  .  ,^(t)  (3.2-2) 

Define  Xj^(l,c,£)  as  follows 


x^(l,C,£)  ■  x^(t,^)  ;  l*l,...,m. 


(3.2-3) 


A  vector  Xj,  whose  1 —  entry  Is  given  by  Xj^(i,t,^),  Is  then  formed 
as  Indicated  below  (See  Fig. 3. 2), 


{xj^d.t,^),  x^(2,t,£),...,  Xj^(n,t,^)}  . 


(3.2-4) 


As  in  previous  techniques  the  number  of  sisnals,  d,  is  estimated. 
Having  an  estimate  d,  a  set  of  (d-fl)  vectors  P  "  l,2,...d+l  Is 
created  by  operating  successively  (p-1)  times  on  X^  with  the  summation 
operatwr,  ^  (*).  For  example,  the  elements  of  Cne  vector  X,  or 

a-3 

>2(3 >£>£)•  J  ■  l,2,...m  are  obtained  as  shown  below: 


X-(j,t,^)  ■  J  X|(l,  C,  ^)  . 

4-j  ^ 


(3.2-5) 
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c,(j)  -  y.  («»-4+i) 

^  a-j 


c,(j)  -  I  c^(a) 


(3.2-6) 


£-j 


Th«  clemencs  of  or  given  by 


x...(j,t,6)  -  I  x.(l,t,Q). 
■"  a-j  “  " 


C3.2-7) 


Using  (3.2-6)  expression  (3.2-7)  becomes 

r  *!,(*•)  (1-1) 

Xjj.l(J.t.0)  -  I  - T  tg(9^)l^^  ^ 


d  s  (t)  a  d  s  (t) 

I  — - — ^  — — 
k-l  [l-g(9^)r  It-l  [l-g(9  ) 


^  (»-j+l)  [g(9^)l' 


k.l 


jIjCjCJ)  I,(9|^))' 


d  ^(t) 
k-l  (l-g(9j^)] 


d-3  I«(VJ' 


I - - -  c  -z^j)  [g(e  )]-. 

k-l  [l-sv9^)]  ^  ^  ^ 
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Therefore, 

b  •  -  in  —  . 

til 

Aaeunlng  b  is  chosen  such  chat  the  Cema  involving  [gCd^)]"  are  negligible 
with  respect  to  the  first  term  in  the  right  side  of  (3.2-8),  x^(J,t,^} 
q  ■  l,2,...,d4-l  and  j  ■  1, 2,..., a  can  be  approxiaated  as  follovsl 


r  i-l 

X  (J.t,9),  I  - lL-_[g(9  )iJ  ^ 

**  k-1  d-gOJ]**  ^ 


(3.2-9) 


As  a  parenthetical  note,  it  is  pointed  out  that  the  above  result  is 
obtained  with  equality  for  a  hypothetical  array  containing  an  infinite 
nuaber  of  sensors.  Then  the  upper  Halt  in  the  suanation  operator  is 
•  and  the  terms  including  (g(6j^>I  do  not  appear.  However,  even  in  the 
case  a  •  «,  this  technique  requires  the  aeasureaents  to  be  weighted  so 
that  terns  involving  the  suaaation  operator  converge . 

Using  the  vectors  1-1,2,...,  d+l,  the  matrix  pencil  M-XN  is  formed 
where  the  rectangular  aatrlces  M  and  N  are  defined  to  be 
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L«t 


E  - 


II  I 

Xj  l2‘"  ^ 

1  1  1 


1  1 

gcep  g(e2> 

g2(e^)  g2(9^) 


g*"^(9,)  g“’^(e,) 


1 1  I 

12  ^3'**  ^+1 

1  1  1 


1 


g(0d) 

g^e^) 


g"-'(9d) 


and 


8^(t)  8j(t)/l-g(9p  ....  s^(t)/[l-g(9^)] 

s^Ct)  S2(t)/l-g(e2) 


d-1 


d-1 


s^(t)  a^(t)/l-g(9^)  ....  8^(t)/[l-g(e^)] 


Using  cha  abova  matrlcaa,  tha  matrix  M  can  ba  dacomposad  as 


d-1 


M  -  EF. 


(3.2-10) 


(3.2-11) 


(3.2-12) 


(3.2-13) 


kl 


Siallarlyt  daflna  th«  dlagoaal  oacrix 


D 


dl«g  [ 


1 

l-g(8j) 


1 

1-8(82) 


1 

1-8(8^) 


]. 


It  follows  thst  N  caa  bo  docooposod  as 

N  "EDF.  (3.2-14) 

Using  tho  docooposlclons  (3.2-13)  and  (3.2-14)  it  is  seen  chat  the 
matrix  pencil  generated  by  this  approach  cakes  the  required  form 
given  In  (2.1-3).  Provided  the  directions  of  arrival  are  distinct,  E 
and  F  are  of  full  rank.  Consequently,  Che  summation  operator  approach 
generates  a  matrix  pencil  Chat  satisfies  the  rsqulresMncs  of  the  pencil 
theorem. 

In  Che  summation  operator  approach  the  11—  entry  d^^  of  the  matrix  D 
Is  given  by 

**11  *  l-g(8^)  • 

In  terms  of  the  generalized  fomulaclon 

“  l-g(e^)  ’  ^  . '*• 

The  values  of  X  which  decrease  the  rank  of  the  pencil  M  -  XN  by  1  arc 
given  by  (2.1-4)  or 

-  U(^^)}"^  -  l-g(8j^)  ;  1-1,2 . d.  (3.2-15) 


Because  the  matrices  M  and  M  are  not  square  In  this  approach,  the 
X^'s  can  be  obtained  by  computing  the  G.E.'s  of  the  square  matrix  pencil 
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glv«a  by  (2.2-1).  Alc«rnatlv«ly.  th«  Gram  macrix  approach  can  be 
uaad  to  obtain  theaa  rank  reducing  numbers. 

Jain's  system  Identification  scheme  has  been  generalized  to  handle 
the  direction  finding  problem.  It  is  shown  that  this  method,  like  ESPRIT, 
la  an  application  of  the  generalized  approach.  This  scheme  makes  use 
of  the  summation  operator.  Being  a  nonsearch  procedure.  It  has  all  the 
advantages  ESPRIT  presents. 

There  Is  still  another  operator,  the  moving  window  operator,  to  be 
discussed  In  greater  detail  In  chapter  four. 


CBAPTER  4 


THE  MOVING  WINDOW  OPERATOR 

Two  oporatora,  cha  phas*  dalay  oparator  and  cha  suomatlon  opara- 
tor,  vara  prasantad  in  chapcar  chraa.  Both  nathoda  vara  fonaulated  In 
tana  of  cha  ganaralizad  approach.  Balng  nonaaarch  procaduraa.  chaaa 
nathoda  poaaaaa  cha  coaputational  and  atoraga  advancagaa  dlacuaaad  aarller. 
Howavar*  cha  ESPRIT  nechod  chat  aakaa  uaa  of  a  phaaa  dalay  oparator  falla 
for  eoharanc  aourcea.  Tha  auoaiacion  oparator,  on  cha  othar  hand,  raqulrad 
a  valghtlng  factor  to  approxlnataly  ellminaca  cerma  chat  would  have  pre¬ 
vented  the  oiacrlx  decompoaiclon  required. 

In  chla  chapter  another  operator  la  preaantad.  Thla  operator,  the 
moving  window,  la  forraulatad  within  tha  ganaralizad  fraoavork  and  la 
ahovn  to  hold  for  coherent  aourcea .  Thla  approach  provldaa  aaynptotl- 
cally  unbiaaad  aatlmacaa  of  cha  anglaa  of  arrival.  Section  4.1  intro- 
ducaa  cha  moving  rectangular  window  for  cha  caaa  of  dacarmlnlsclc  alg- 
nala.  In  aubaactlon  4.1-1  cha  caaa  of  a  alngular  algnal  covariance 
matrix  la  analyzed.  Tha  rectangular  moving  window  la  chan  applied  to 
cha  zero-mean  random  algnal  caaa  In  aubaactlon  4.1-2.  In  aecclon  4.2 
cha  moving  rectangular  window  la  ganaralizad  to  a  rectangular  planar 
array.  Sacclona  4.3  and  4.4  praaanc  two  wall  known  machoda,  Prony'a 
and  Plaaranko'a  machoda  [42-46].  .  Thalr  ralatlcnahlp  to  tha  moving  rec¬ 
tangular  window  la  developed.  A  comparative  performance  baaed  on  com¬ 
puter  almulatlon  of  cha  chraa  machoda  la  given  In  aecclon  4.5. 

4.1  Tha  Moving  Bactangular  Window 

In  this  aactlon  a  nonaaarch  procedure  la  praaancad  for  azlmuch-only 
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DOA  •selmatlon  of  far  flald  point  sourcas.  Consldar  a  llnaar  array 
eooposad  of  a  Idanclcal  oanldiractlonal  sanaors  with  uniform  sensor 
spacing  0.  Aasuaa  chare  are  d  ^  a/2  narrovbaud  deterministic 
sources  located  at  azlauchal  angles  3^,  1  ■  1,2,. ..fd,  which  are  Im' 
pinging  on  the  array  as  planar  wavefronts  and  emitting  signals  whose 

coaplax  envelopes  are  denoted  by  s^(c),  ]c>l,2 . d.  The  signal 

received  at  the  1^  sensor  Is  the  superposition  of  d  Impinging  wavefronts 
plus  zero-mean  additive  noise.  Define  relative  response 

of  the  1~  sensor  to  the  source,  oIq  to  be  Che  center  frequency  of  each 
of  the  spatial  sources,  c  to  be  the  speed  of  propagation  of  the  plane 
waves,  and  n^(t)  to  be  the  additive  noise  of  the  1~  sensor.  Then  the 
signal  received  at  the  1—  sensor  can  be  expressed  as 


d 

y  s.(t)  a/0.)  >»'n.(t);  l»l,2,...,n. 
k«l  “  ^  ^ 


Let 

<|)j^  -  ~  D  sin  ;  k«l,2,...,d. 


(4.1-1) 


Vftien  using  a  linear  array  of  omnidirectional  sensors,  the  relative  res¬ 
ponse  (senser  1  being  the  reference  sensor)  of  the  1~  sensor  to  the 
source  Is  given  by 

(i-l) 

•i^V  "  *  •  (4.1-2) 


y 


1 


(t,^)  can  Chen  be  rewritten  as 


d 

y.(c,e)  -  ^  j,  (t)  exp{{j<^  }  (i-l)]+n  (t) 
1  k-l  “  “  ^ 


•  I  s.(t)  exp[(j~D  sin  9.)  (1-1)]  +  n/t). 
kii  k  c  k  1 


(4.1-3) 
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Th«  axp^ctsd  valu«  of  yields  x^(t,^)  given  by 


d 

-  y  i“l-t2, . . . ,m. 

k-1  ^  “ 


Define  Che  rectangular  window  co  be 

a^,(n) 


tnguLar  winaow  co  i 

fl  1  1  a  ^  N 

1^0  otherwise. 


Given  Che  ouaber  of  sources  d  and  Che  m  averaged  data  points  x 
a  sec  of  (<tfl)  vectors  ^  is  created  where  the  conponencs  of  Z 
Che  (n-d)  values  of  che  sequence 

Xj^(e,^)R^  (l-n”*-!);  a  ■  l,2,...,d  +l 

n£i£«  +  n-  d“l. 

In  particular,  (see  Fig.  4.1) 

-  {x^(t,£),  XjCt,^),...,  x^  j(t.9)} 

X^  •  {x2<t,£),  Xj(t,£),..., 


^+1* 


The  matrix  pencil  H-XM  is  then  formed  where 


(4.1-4) 


(4.1-5) 


i(t.9), 

arc 


(4.1-6) 


(4.1-7) 


(4.1-8) 
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D«fin«  th«  aatrlcas 


A  » 

n 


•n<®l> 

‘n<®2> 

...  a  (9 
n 

*n.l<®2> 

•  • •  *n+l 

(4.1-9) 


E  •  A 


-  A, 


0*1 


(4.1-10) 


and 


♦  -  <il*8  d22 . 

J02  ^^d 

■  dlag  [a  ,  •  “] 


(4.1-11) 


Also,  dsfln*  ch«  vector 


£  "  *2 . *d^‘ 


It  can  be  shown  that  the  vector  X  can  be  expressed  as 


X  >  A  S;  n*l,2,...,d'fl. 
n  •“ 


(4.1-12) 


Observe  fron  (4.1-2)  that 

/__n  j0^(n-l) 

-  {«2(V^^  -  {e  ^  }. 


It  follows  from  the  above  observation  that 
A  -  A.0 

Q  I 
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Consequently , 


Def ine 


X  -  A  S  -  E*  S. 

—no—  — 


j<l>.  j(d-l)(t. 

Sj^e  ....  s 


2  .1(d-l)«2 


j4> 

a^e  ....  82 

*  J<*>d 

S  ,6  *  •  •  •  S 

a  a 


j(d-l)^. 


The  matrix  M  can  then  be  expressed  as 


(4.1-13) 


(4.1-14) 


M  -  EF 


(4.1-15) 


where  E  and  F  are  given  by  (4.1-10)  and  (4.1-14),  respectively. 
In  a  similar  manner  N  can  be  expressed  as 


N  -  E'SF. 


(4.1-16) 


Using  (4.1-15)  and  (4. 1-16), the  matrix  pencil  M-\N  can  be  decom¬ 
posed  as 

M  -  XN  -  E(I  -  X$)F.  (4.1-17) 


This  decomposition  satisfies  the  requirements  of  the  pencil 
theorem.  As  long  as  the  directions  of  arrival  of  the  signals  are 
distinct  and  the  separation  is  less  than  half  a  wavelength,  the 
columns  of  E  are  linearly  independent.  Thus,  the  matrix  E  is  of 
rank  d  as  required.  By  inspection,  the  matrix  F  is  also  of  rank  d. 
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H«ac«.  ch«  values  of  X  for  which  the  rank  of  tha  matrix  pencil  M-XN 
dacraaaas  by  1  art  given  by 


-I 

*  ‘^kk  "  *  •  *  1,2,.. .,d. 


(4.1-18) 


The  values  of  X  can  be  cosiputed  by  means  of  a  Gramnian  approach. 

The  rank  reducing  X's  are  the  values  of  X  for  which  the  columns  of 

the  matrix  pencil  M-XM  become  dependent.  Define  the  inner  product  of 

the  vectors  ^  and  ^  by 
— i 


(4.1-19) 


where  (*)°  denotes  the  Hermitian  operator.  To  check  the  dependence  of 
the  set 


<il  -  *2  -  % . id  -  'W- 


(4.1-20) 


construct  the  Granmian  matrix 


^•2™  XXj  ,X^— XX^^  ^X^— XXj  *  ^5.3^  ...  ^»2  *  XX^  d  *  ^^d+1^ 


(4.1-21) 


A  polynomial  P(X)  whose  xeros  are  the  rank  reducing  numbers  can 
be  obtained  from  the  equation 


det  (G)  •  0. 


(4.1-22) 
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An  nlcarTintlvn  approach  to  cha  cogq>ucaClon  of  the  X  is  glvan 
in  sacclon  2.2.  It  is  shown  that  tha  rank  reducing  X's  of  H-XN  are 
also  the  generalized  eigenvalues  of  the  matrix  pencil  given  by 

-  XM^'n.  (4.1-23) 


Again,  once  the  G.E.'s  of  (4.1-23)  are  known,  (2.1-2)  can  be  used 
to  compute  the  DOA's. 

In  suonary,  given  M  snapshots  and  the  nuaiber  of  sources  d,  the 
following  algorithm  is  proposed.  It  is  assumed  that  the  N  snap¬ 
shots  are  taken  at  a  rate  sufficiently  fast  such  that  s^(c),  k>l,2,...,d, 
remain  approximately  constant  over  the  M  snapshots. 


ALGORITHM  1 

Step  1;  Given  N  snapshots,  form  the  averaged  data 

1  ” 

*i^t.i)-^  I  y  (t  ,6);  i-l,2 . . 

”  j-1  ^  J 


(4.1-24) 


Step  2;  Form  the  vectors  i«l,2,...,  d+l,  in  (4.1-7)  by  approximating 


x^(t,^)  with  x^(t,£). 


Step  3;  Form  the  matrices  M  and  N  whose  columns  are  given  by  ^  and 
,  i  ■  1,2 . d,  respectively. 

— l-rl 

Step  4;  Form  the  matrices  m”m  and  M^N. 

Step  5;  Compute  the  Generalized  eigenvalues  of 

-  Xl^. 
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Th«M  G.E.'s  ar«  known  to  b«  of  the  form 

0). 

\  -  sxp  {-j  —  D  8ln  9j^}. 

St«p  6;  Find  ch«  DOA's  using  ths  rslatlon 

“  sresin  {j  *  l»2....td. 


(4.1-25) 


(4.1-26) 
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4.1>1.  Si 


ar  Signal  Covarlanea  Matrix  C«s«: 


An  ioportanc  proparty  of  cho  aovlng  window  approach  Is  chat 
Cha  propoaad  algorithm  doas  not  fail  If  cha  sourca  signal  covariance 
matrix  becomes  singular.  Singularity  of  the  covariance  matrix  may 
be  due  to  sturt  Jammers  using  coherent  signals t  the  presence 
of  multipath  phenomenon,  etc.  To  verify  this  property  consider, 
once  again,  the  matrix  pencil 


M-XM  -  E(1-X*)F 


where  E,  and  F  are  given  by  (4.1-10),  (4.1-11)  and  (4.1-14), 
respectively.  Mote  chat  the  matrix  E  Is  of  rank  d  as  long  as  the 
directions  of  arrival  ar.e  distinct  and  the  sensor  separation  Is  less 
than  X/2.  The  matrix  F  Is  also  of  rank  d  even  In  the  presence  of 
coherent  sources.  The  diagonal  matrix  (I-X4),  whose  11^  entry  Is 
(1-Xe  ),  Is  of  rank  d  provided  A  f(  e  ;  l«l,2,...,d.  Therefore 

the  requirements  of  the  pencil  theorem  are  satisfied  even  when  the 
signals  are  coherent.  Because  d  rank  reducing  numbers  exist,  even  In 
Che  case  of  coherent  sources.  It  follows  chat  the  moving  window 
approach  Is  applicable  even  when  the  signal  covariance  matrix  is 
singular . 


4.1-2  Zero-Ifcan  Itodom  Signals  Case; 

In  Che  moving  window  approach  presented  In  section  II  the  sources 
were  assumed  to  be  deterministic.  Had  they  been  random  with  zero 
mean,  x^(t,_9);  i«l,2,...,m  would  have  been  zero.  The  matrices  M  and 
N  In  (4.1-8)  would  Chen  be  Identlally  zero  and  would  be  of  no  use  In 
determining  the  OOA's.  In  this  section  it  is  shown  that  a  modification 
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of  ch«  oovlag  trladow  approach  is  applicable  whaa  cha  aoa-zsro  amma 
sisnal  asauapclon  does  not  hold. 


Aaauoa  d  random  narrowband  stationary  sources  with  zero  mean  and 
canter  frequency  uJq.  Given  the  vectors  1  •  d+l,  where 

y2(t.6)f-t  y^^^j(t,9_)} 

I2  “  ty2(t.6).  73  (t  .9) . 

The  Inner  produce  <  •,•  >  is  now  defined  to  be 

<4.  lj>  -  EtYj-Y^l.  (4.1-27) 

Define  the  matrices 

(4.1-28) 

(4.1-29) 


<l2* 

<Yj,  Y^>  <Y3, 


•  V 

.  <ly  Y^> 


- 


* 

<Y, ,  Y->  ... 

<Y,  , 

—1 

—1 

—1  —2 

—1 

<^1. 

<Y,,  Y,> 

<Y-, 

—2 

—1 

—2  —2 

—2 

V 

<Y^,  Y2>  ... 

"2d* 
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A  - 


1 

1 

1 

a 

j<>2 

6  a  a  a 

a 

a 

* 

32*2 

6  •  •  a 

a 

m, 

a 

a 

a 

(n-d-l)(tj^ 

a 

j(n-d-l)(>2 

a 

j(B-d 

a 

(4.1-30) 


(4.1-31) 


and 

«  ■  dlag  [a  ,  a  , ....  a  ] , 
Also,  daflna  cha  vectors 

T 

^  *2»***» 


^  “i+i»**’'\+(m-d-l)^ 


It  follows  that  can  ba  axprassad  as 

Y^  -  A* 

Using  (4.1-32),  tha  inner  product  <^,  bacoaes 
<Y.  ,  !,>  -  <Ai  +  N.  .  S  +  ^> 

-  E((At^^"^^  S  +  ^)®(A*^’‘“^^S  ^)1  ■ 


Aasuiaa  tha  zero  aaan  noise  conponants  with  variance  a  to  ba  uncorrela- 
tad  from  sensor  to  sensor  and  independent  of  cha  signals.  Then 


(4.1-32) 
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‘’‘•''si  +  (-AXJ^ 


I  ^  k. 

i  -  k.  (4.1-33) 


Uit  I,  b«  d«£ln«d  as 

0  1  •  .  .  •  0 

0  0  1  ...  0 

.  .  • 

. 1 

. 0 

Assualr.g  chs  oolss  powsr  to  bs  kao%m,  the  matrix  pencil  M-XN  is 
defined  as 

M-XN  •  (M^  -  (»-d)a^I)  -  X(Nj^  -  (n-d)o^ip. 


Observe  that  the  kf/^  element  of  M  is 

EIS“  a“a«<‘-*>S) 


(4.1-34) 


whereas  the  kl'^  element  of  N  is 

S].  (4.1-35) 

Let  us  now  get  a  closed  form  expression  for  the  expectation  given 
in  (4.1-34).  To  obtain  this  result  we  perform  the  matrix  multiplica¬ 
tions  involved  in  (4.1-34). 
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.-d-l  J(*  ,“♦,)! 
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jU-lH,  ■-d-l  J(k-l)<|>_  m-d-l 


sV("'i)  a“a  s  - 


j(k-4)«,  j(k-l)^.  n^-l  J(<t,-<l>,)i 

*8,(B-d)  •  ^  +8  *s,  •  ^  (1  +  I  e  ^  ^  ) 

j(k-l)*.  »-d-l  j(0. -<(»,) 

+...  +  8  *8.  ^8  ^  (1  +  I  e  M  + 

^  ^  i-1 

-ja-l)<D2  J(k-1)<D  .-d-l  Jf.  .w  j(k-£)^ 

*8^  8  ^8  Ml  +  I  8  +  82*82(m-d)8  ^ 


8^*8^ (nrd)  8 


■►. .  .+■  »2**d  * 


-J(k-l)d,  j(ll-l)<l>.  Il-d-1 

^8  ^  (1  +  y  8  M  + 

i-1 


83*818 


-j (4-1)0,  J(k-1)0,  »^-l  j(0,^,)i  *  -j (4-1)0,  j(k-l)0, 

^  Ml  +  I  8^3+  8:8,  8  M  ^ 
1-1 


m-d-1  J(0,-0i)i 

(1  +  I  8  2 

i-1 


)+...+  838^  8 


*  -j  (4-1)0,  J(k-1)0^  «-d-l 

,8.  8  8  (I  +  y 

J  °  i-1 


*  *d  *1  • 


-j(4-l)0.  j(k-l)0,  m-d-l  j(0,-0.)i  *  -j(4-l)0.  j(k-l)0, 

8  Ml  +  I  8  ^  **  )+  8^8  8  M  2 

i-1 


m-d-l  j:0,  -0.)1  *  j(k-Q0 

(1  +  y  8  )+...+•  8  8.(m-d)  8 


Th8  8bov8  8qu8Clon  can  ba  vritcan  in  a  much  more  compact  form  using 
cha  following  auatltutiona .  Lac 


F  •  1  +  I 


"^-1  nK 


p  q 


SH,H(4.1)  ,(k-l)g 


d  d 


-  y  y  8*  8 

q.l  pil  * 


-J  (4-1)0.  j(k-l)0„  tD-d-1  j(0_-0  )i 


P(l+  I  8  P  ) 
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To  slopllfy  further  d«£ln*  th*  matrices 


and 


U  - 


1 

1 

1 

iOl 

J«2 

J^d 

e 

e 

e 

J2d 

e  ^ 

e 

• 

• 

• 

• 

• 

j(«l-l)^2 

j(d- 

• 

e  ^ 

e 

1“ 

^l^l 

^12^12  •  •  • 

^d" 

^21^21 

*12^22  •  •  • 

S2d^ 

• 

• 

• 

• 

s 

• 

« 

• 

• 

^dl^’dl 

^d2^d2  *  •  • 

Sdd". 

The  matrix  M  can  than  be  decomposed  as 

M  -  OVU^  . 

Similarly,  the  matrix  H  can  be  decomposed  as 

N  -  OVfV. 


(4 


(4, 


(4, 


Thus,  as  required  by  the  pencil  theorem,  the  matrix  decomposition  of 
M-XN  is  given  by 

M-XH  •  UVU®  -XUV«®U® 

-  OV  (I  -X*®)0®. 


1-37) 


1-38) 


1-39) 


1-40) 
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The  matrices  UV  and  U«  are  of  rank  d  as  long  as  Che  directions  of  arrival 
of  Che  signals  are  dlsclncc.  Hence,  che  values  of  X  for  lAlch  Che  rank 
of  Che  nacrlz  pencil  M-XN  decreases  by  1  are  given  by 

JOi 

X^*e  ^  1*  l.l.***.  d« 

Note  Chat  che  matrices  UV  and  are  of  rank  d  even  when  che  sources  are 
coherent.  Therefore,  as  in  che  nonzero-^ean  signal  case,  che  approach  does 
noc  fail  in  che  presence  of  coherent  sources. 

The  above  choices  for  M  and  H  are  noc  unique.  To  demonstrate  che 
flexibility  of  che  pencil  cheorem,  a  second  choice  is  now  considered. 


(4.1-41) 


1  k  -  i 


0  otherwise . 


Using  (4.1-32),  and  become 

+  AifS  A®  +... 

...^A<f  W-i)s  + 

Nj^  -  EfA*  S  A®  + 

...  A^<‘^>isV<J-^>A«Wa2l[l-6,  .1. 

j » ^ 


(4.1-43) 


(4.1-44) 


(4.1-45) 
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Lac 

M  -Mj  -  ahll  - 

and 

N  - 

Daflna  cha  nacrlces 

S  -  E[S  , 

T  -  [S  $S  +. ..+ 
-  [I  +  ♦  +...+ 


(4,1-46) 

(4.1-47) 

(4.1-48) 

(4.1-49) 


and 

2  -  A®.  (4.1-50) 

Iha  macrlx  H  can  chan  be  decomposed  as 

M  -  ATZ. 

Similarly,  N  can  be  decomposed  as 
N  -  A  ♦!  2. 

Consaquencly,  as  requirad  by  Che  pencil  Cheorem,  che  macrlx 
pencil  M-/.N  can  be  decomposed  as 

M  -  XN  -  A(I-Xt)T2  .  (4.1-53) 

However,  Chls  approach  falls  whan  che  sourcas  are  coherenc  since 
cha  macrlx  T2  Is  reduced  In  rank  for  chls  case. 


(4.1-51) 


(4.1-52) 
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4.1-3  Pr«flltTing 

Th«  wing  window  approach  yields  correct  results  la  absence 
o£  noise  assunlng  sufficient  ntmerlcal  accuracy.  Therefore*  it  la 
reasonable  to  expect  that*  If  the  SNR.  In  the  data  can  be  Improved 
by  prefiltering*  then  the  performance  of  our  approach  would  Improve. 

The  question  is*  "Ooca  prefilterlng  prevent  application  of  the  matrix 
pencil  approach?"  In  this  section  we  apply  a  frequency  selective  fil¬ 
ter  at  each  channel  (see  Fig.  4.2)  and  show  that  the  matrix  pencil 
decomposition  Is  not  disturbed. 

Assume  that  a  filter  with  Impulse  response  h(c)  Is  applied  to 
each  channel.  The  filtered  output  of  the  1~  channel  Is  given  by 

•  I  h(t-T)dT  (4.1-54) 

o 


where 

y^C". 


[o*T]  la  the  observation  Interval. 
6)*  u. (t*  9)  becomes 

-  1  j 

0 


d 


Using  the  expression  for 


u^(t*6) 


\  s  (■0a^(9j^)h(t-T)  dx 
k*l 


rT 


+ 


4) 


n^(T)h(t-x)dT 


d  T 

-  y  s.(T)h(t-T)  dx 

kil  ^ 


t  T 


n^(t)h(t-x)dx  . 


(4.1-55) 
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L«t 


*1^(0 


8^(T)h(t-T)dT 


(4.1-56) 


and 


n^  <c)  -  I  n^(T)h(c-T)dT  . 
1 


(4.1-57) 


Tha  axprasaion  for  u^(ttd_)  can  than  ba  rawrlttan  In  carma  of  z^vt) 


and  n  (c)  aa 
“t 


u.(t,i)  •  I  2it)a  (6  )  n  (t). 
1  k-l  “  ^  *  “l 


(4.1-58) 


Asaualng  tha  oolaa  to  ba  zaro  aaan  and  tha  aignala  to  ba  dater- 
■InlaclCt  tha  axpactad  valua  of  u^(t,£)  la  glvan  by 

■  «1  %<'>•!< V  *  V<'» 


k«l 


I  *k(t)a,(9.)  +  Etn  (t)l  . 
k-1  *  z 


(4.1-59) 


By  uaing  a  ractangular  window,  a  aat  of  (d+1)  vactora  la 
craatad.  In  particular. 

vj  -  {v^(t,£),  V2(t,0.) . 


T 

V 

^+1 


.,  v_(t,6)}. 

m 


(4.1.60) 
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Th«  aatrlx  p«iicll  M-XM  it  formtd  lAtrt 


M  • 


I  I  i 

li  ^2  •  •  •  ^ 

1  1  1 


M  - 


[nr 

V  7  7 
^  -3 

1  1  1 


L«C 


2^  -  [*j(t) . *d(t)l 

Tht  vtccor  can  b«  txprttttd  tt 


whtrt.  at  in  (4.1- 

13).  (4.1-10),  (4.1-11), 

Otflna 

A,  -  £♦ 

cht  matrix 

r 

z^(t)  *j(t)t 

j(d-l)0,  “ 

*j(t)  C 

J^2 

J(d-1)(> 

*2(t)  e  ^ 

« 

V 

F  - 

• 

1 

• 

• 

• 

*^(c) 

*^(t)t  ®  .  .  . 

j(d-l)d. 

Zj(t)  e 

Tht  ntcrlcts  M  and  N  can  bt  txprtsscd  as 
M  •  EF  and  N  -  E*F  . 


(4.1-61) 


(4.1-62) 


(4.1-63) 


(4.1.64) 


(4.1-65) 


71 


lhu««  tb«  aacrix  pencil  M-XN  can  ba  daconposad  as  requlrad  by  tha 
panel!  thaorui.  Spaclflcallyi 

M  -  XN  -  EF  -  Xe^F 

-  E  (I  -  X«t)F  .  (4.1-66) 


Tha  aatrlz  dacoaposltlon  ls»  charafora*  prasarvad  without  altering 
tha  diagonal  matrix  4  idian  nultlchannal  filtering  Is  added  to  the 
moving  window  approach.  As  before •  the  rank  reducing  numbers  of  tha 
matrix  pencil  M-Xn  will  still  contain  tha  necessary  Information  to 
astlmata  tha  OOA's. 

Tha  purpose  for  using  a  frequency  selective  filter  In  each  channel 
Is  to  Improve  the  SNR  In  the  entries  of  the  matrices  M  and  N.  A  simple 
example  Is  discussed  next  to  Illustrate  the  procedure.  The  power  spectral 
density  P(u)  of  the  noise  Is  assumed  to  be  constant  up  to  some  high  cut¬ 
off  frequency  such  that 


P(u) 


No/ 2 


— u)  <  U>  <  U) 

c  c 


|(i)| 


(4.1-67) 


Assume  the  desired  signal  with  the  largest  bandwidth  to  have  bandwidth 

AB.  Let  P  denote  the  average  power  of  the  k^^  desired  signal.  Suppose 
*k 

now  that  the  received  signals  are  put  thru  a  rectangular  bandpass  filter 
with  center  frequency  and  bandwidth  AB.  The  t.lgnal-coHnolse  ratio  at  the 
Input  of  the  filter  Is  given  by 

d 

J.  \ 

o  c 
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wfall*  ch«  •Igaal  to  oolao  ratio  at  tha  output  of  tha  flltar  la  glvan  by 


d 

o 

Daflna  tha  algnal-to-oolaa  ratio  laprovaaant  factor  to  ba 

(SMI) 

^•Ts5r57  • 


For  tha  aoving  wlndov  oparator  and  tha  signals  posad  In  this  asasipla, 
tha  SMt.  of  tha  antrlas  In  tha  aatrlcas  M  and  N  ara  Incraaaad  by  tha 
factor 


I 


(4.1-71) 


Ona  can.  tharaforat  loprova  tha  parfomanca  of  tha  aovlng  window 
approach  by  praflltarlng  tha  signals.  Tha  flltars  Introducad  into 
aach  channal  do  not  disturb  tha  oatrlx  daconpositloo  of  M-A.N  as  required 
by  the  pencil  thaoram.  Tha  SMI  is  Increased  In  the  entries  of  M  and  N 
and  ona  can  expect  aora  accurate  estlaatas  of  the  DOA's. 
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Thtt  sch«M  propo««d  in  snctlon  4.1  also  lands  Itsalf  co  a  rectangular 


planar  array  of  sensors  with  uniform  spaclngs  d  >  d  In  the  x>  Y  dlrec- 

X  y 

cions  respectively  (See  Fig.  4.3). 

The  angular  locations  of  d  sources  are  again  to  be  estimated.  These 
estimates  are  to  be  obtained  from  measurements  collected  at  the  sensors 
of  a  rectangular  planar  array.  Assume  this  rectangular  planar  array 
to  have  p  sensors  in  the  X  direction  and  q  sensors  in  the  Y  direction. 
Observe  chat  one  can  chink  of  this  array  as  being  comprised  of  either 
p  columns  of  length  q  or  q  rows  of  length  p.  In  this  section  ve  view 
the  array  as  being  composed  of  q  rows  of  length  p.  To  be  able  to 
locate  d  sources  for  this  case,  q  has  to  be  greater  chan  or  equal  to 
(dH-l)  while  p  must  be  greater  chan  or  equal  to  d.  Define  i  and  1 

y 

to  be  Che  unit  vectors  along  the  Z  and  Y  axes.  The  position  of  the 
im~'  sensors  is  given  by 


Uc  -k  be 
— r 


The  signal 

on  — 

where  Z  • 
-mn 


^  -  (m-l)d  1  +  (n-l)d  i  . 

X-X  y-y 


(4.2-1) 


Che  direction  of  propagation  for  the  signal  where 


k  •  cos  9  1^  +  sin  9„  1  . 

-r  r  -X  r  y 


ch 


(4.2-2) 


received  at  the  nai  sensor  can  therefore  be  expressed  as 
d 


s^(t)  exptj  +  n^(t); 


(4.2-3) 


k  denotes  the  inner  product  of  the  two  vectors  Z  and  k  , 
— r  -on  -r 
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SOURCE 


-k 


Fig.  4.3  Rectangular  planar  array 


and  a^(t)  danotas  cha  addltlva  noisa  at  the  on —  sanaor.  Using 
(4.2-1)  and  (4.2-2),  cha  tnnar  product  ^  results  in 

Z •  k  •  [(n-l)d  1  +  (n-l)d  1  l.[Cos  9_  1  +  sin  9_  i  1 

-Tsn  m  y~y  r  r  -7 

■  (at-l)d^  cos  9^  +  (n-l)dy  sin  0^.  (4.2-4) 

Insarcing  (4.2-4)  Into  the  expression  for  y  (c,6),  we  obtain 

m  • 

y^(C,^)  »  S^(C)  exp  [j~  (a-l)d^  cos  9^]  exp  [J  ^  (n-Dd^  sin  9^.1 

+n  (t).  (4.2-5) 

nn 

For  example,  the  signals  at  the  sensors  In  the  first  rov  of  the  array 
are  given  by 

d 

Yll  •  I  s  (t)  +  n  (t) 
r-l  ^  “ 

d 

721  ■  I  X  '''  ”21^^^ 

•  d  - 

ypl  •  I  s  (t)  exp  [j  ^  (p-l)  d^  cos  9j.l  +  n  (t) . 
r«l 

Assuming  the  signals  to  be  deterministic  and  the  noise  to  be  zero  mean, 

Chen 


X  (t,e)  “Ety  (t,0)l 
on  •  nm 

■  ^  s  (t)  exp  tj  d  (m-1)  cos  9  ]  exp[j  d  (n-l)sln  6  ].  (4.2-6) 

IP  A  ^  • 
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Oefliw  Che  mecrlces 


77 


and 


whara  E  and  F  ?ra  glvan  by  (4.2-U)  and  (4.2>1S),  raspectlvaly.  In  a 
sinllar  mannar,  N  can  ba  axprassad  as 

N  ■  E'tF  .  (4.2-17) 
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Using  (4.2-16)  and  (4.2-17)*  chs  astrlx  psncll  M-Xtt  can  ba  dcconpoaad 


M  -  XN  -  EF  -  Xe6F 


-  E(1  -  X4)F. 


(4.2-18) 


This  dacoapoalclon  satisfies  cha  raquiranants  of  Cha  panel!  chaoreo. 
Since ,  in  general ,  M  and  N  are  rectangular  matrices,  either  one  of  the 
methods  suggested  in  section  2.1  can  be  used  to  compute  the  generalized 
eigenvalues.  Observe  that  q  generalized  eigenvalues  exist  in  this 
case  but  only  d  of  them  will  lie  on  the  unit  circle.  These  d  G.E.'s  are  of 
the  form 


^i  •  ‘‘ii  ’  •  »  . 


(4.2-19) 


The  remaining  (q-d)  G.E.'s  are  at  the  origin.  Thus,  once  the  d 
nonzero  G.E.'s  are  known,  the  DOA.'s  can  be  calculated  from 


-  arcsin  {j  —  In  '{'2^^  »  i“l»2,...,d.  (4.2-20) 

y 
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4.3  Prony**  Method 

In  this  ssctlon  Prony's  oschod  [44]  Is  prsssncsd  in  the  coatext 
of  array  signal  procasslng.  Wa  chan  show  that  Prony's  algorithm 
Is  ralatad  to  tha  moving  ractangular  window  prasantad  in  sactlon  4.1. 

Prony's  machod  Is  usad  for  nodallng  data  of  aqually  spacad  samplat 
by  a  linaar  comblnacion  of  axponantlals .  Tha  original  procadura  by 
Baron  da  Prony  flttad  axactly  d  axponantlals  to  2d  data  maasuramants. 
For  tha  casa  in  which  only  anapproxlmata  fit  of  d  axponantlals  to  a  data 
sat  of  a  saa^las  Is  daslrad.  wfaara  m  >  2d«  a  Isaac  aquaras  astlaatlon 
procadura  la  usad.  Tb.ls  procadura  Is  callad  tha  axtandad  Prony  machod. 
Tha  modal 


d  ‘'k(i-l) 

y.(t)  •  I 

^  k-l  * 


+  n^^Ct);  1-1, 2,..., m 


(4.3-1) 


Is  to  ba  usad  In  tha  axtandad  Prony  mathod  for  approximating  tha 
maasurad  data  y, (t) ,  y-(c),...,y  (c).  Tha  problam  is  to  astlmaca  tha 
paramatars  Z^;  k-l,2,...,d.  This  formulation  Is  similar  to  chat 
usad  in  array  signal  procasslng  whare  tha  problam  is  to  astlmaca  the 
angular  locations  of  d  sourcas,  a^(t}  reprasancs  tha  complax  envelopes  of 
the  signals,  and  are  purely  imaginary  parameters  (the  damping  factors 
are  equal  to  zero) .  ^  given  by 


^k  “  ^’*’k  “  ^  c  ° 


(4.3-2 


where  ui  ,  c,  0  and  6.  are  defined  in  sactlon  4.1.  J  is  tha  square  root 
o  ^ 

of  minus  one.  Because  of  tha  dapandanca  of  the  modal  on  the  angles  9^. 
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(4.3-3) 


k"l,2 . da  ch*  aMsuraMnts  y^(t)  ar*  written  as  y^(t,^).  Thus,  the 

dual  to  equation  (4.3-1)  Is 


yt<t.e)  -  I 
^  k-l 


•k(t) 


k(l-l) 


+  n^(t) 


which  Is  Identical  to  (4.1-3).  Aasualng  the  signals  to  be  deterministic 
and  the  noise  to  be  zero  mean,  the  expected  value  of  y^(c,£)  yields  x^(t,;9) 
given  by 

X  (t,^)  -  I  a.(t)  e  ;  1-1,2 . m.  (4.3-4) 

*■  k-l  “ 


Finding  the  values  k>l,2,...,d  that  minimize  the  squared  error 

la  a  difficult  nonlinear  least  squares  problem.  An  alternative  suboptlsaim 
solution  was  suggested  by  Prony. 

The  key  to  the  Prony  technique  Is  to  recognize  that  (4.3-4)  Is  the 
homogeneous  solution  to  a  constant  coefficient  linear  difference  equation, 
the  form  of  which  la  developed  next.  Expanding  the  expression  for  x^(t,^) 
we  obtain 

x^(t,^)  -  s^(t)  exp[j(>j]  (1-1)  +  3^(1)  exp[j(fr2l  (i-U 

+. ..+  Sj(t)  exp(j0^] (1-1) .  (4.3-5) 


Using  (4.3-2),  (4.3-S)  becomes 


x^(t,9)  ■  Sj(t)  H 


...+  s^(t)  Z^ 


(1-1) 


(4.3-6) 


Evaluating  (4.3-6)  at  l>l,2,...,m,  we  obtain 
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■•“  •2<t)Z2+*-*+ 

(4.3-7) 

x^(c,e)  +...+  g^(t) 

Aa«ua«  ctuit  ch«  coapl«x  «zpao«aclalg  Z^.  k-lt2....,d  of  (4.3-2)  «r«  ch« 
roocs  of  th«  algabralc  •quaclon  (also  known  as  cha  pradlction-arror  filter 
polynoolal)  given  by 

+  S2Z  +  8^2^  +•••+  -  0  (4.3-8) 

where  la  arbitrarily  set  equal  to  1.  In  order  to  detezvlne  the 
coefficients  Bp  S2>«**>Sj>  the  first  equation  of  (4.3-7)  is  nultlplled 
by  Bp  the  second  one  Is  nultlplled  by  B2>  and  the  one  Is  nultlplled 
by  B^  and  finally  the  (d-fl)tb  Is  multiplied  by  B^^^  •  I.  This  results  In 
the  set  of  equations 

SjX^(t,^  ■  +...+  Sia^(t) 

82*2^^’-^  ■  B2Sj(t)Zj+B2S2(t)Z2  +...+  B2S^(t)Z^ 

■>•  8d*2‘'>*2'‘‘‘  *•••*  ®d*d<'>*d'‘ 

*d+l^^*^  -  s^(t)Z^  +  S2(t)Z2  +...+  s^(t)Z^  . 
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Adding  th«  above  aquaclons  and  using  cha  face  chat  Z, ,  Z. are  chc 

i  Z  a 

roots  of  (4.3-8).  «a  have 

(c.e)  - 

s^(c)  [Sj  +  SjZ^  +...  +8^  zj"^+  zj] 

+  s,(e)  tS.  +  8,Z,  +  ...  +  8.Z^“^  ,d, 

2  -  22  d  2  +  Z“] 

+  2 

« 

« 

+  a^(t)(Sj  +  822^  +  ...  +  Sj  2^“^  +  Z^  ]  -  0. 


Similar I7,  a  sat  of  (m-d-l)  additional  aquations  having  the  same  form  is 

obtained  by  successively  starting  with  the  second,  the  third . the  (m-d)ch 

equation.  For  convenience,  ve  vrite  x^(c,^}  as  in  the  following  matrix 

equation.  The  set  of  equations  obtained  following  the  procedure  is  written 
in  matrix  form  as 
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(4.3-10) f  which  coaprlsAS  (a-d)  •qu«tlona  In  d  unknowns,  can  bs  solved 


by  using  s  Issst  squares  approach.  Having  found  the  S^'s,  the  algebraic 
equation  (4.3-8)  is  forasd  and  its  zeros  cooputed.  Once  these  zeros, 

Z^:  k^l,2 . d,  are  known,  expression  (4.3-2)  is  used  to  obtain  the 

directions  of  arrival  1a*l,2,. . .  ,d. 

The  relationship  between  Prony's  aethod  and  the  aoving  rectangular 
window  technique  is  deaonstrated.  These  aethoda  are  related  in  the 
sense  that  the  inforaation  needed  to  estiaate  the  angular  positions  of  the 
d  sources  is  obtained  by  exaaining  the  dependence /independence  of  a  set 
of  vectors  formed  froa  the  averaged  vectors  i*l,2 . drt*l  given  in 

“X 

(4.1-7). 

Sacall  that  the  aoving  window  approach  exaaines  the  dependence  of  the 
set  of  vectors  i*l,2,...,d  given  in  (4.1-20).  On  the  other  hand, 

the  algebraic  equation  (4.3-8),  whose  zeros  are  given  by  Z^,  k*l  ,2,...,d  in 
(4.3-2),  can  be  derived  by  exaaining  the  dependence  of  the  vectors 
:  i"l,2,. . .  ,d4-l,  such  that 

^252  2^+1  "  -•  (4.3-11) 

In  view  of  (4.1-10),  (4.3-11)  can  be  %rritten  as 

SjES  +  S2E*S  +  ...  +  i  ■  0  (4.3-12) 

or,  equivalently, 

EtS^  82'^  +  ...  +  S  -  0 .  (4.3.13) 
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For  cquAllty  to  hold  In  (4. 3*13).  th«  matrix  equation 


6i  +  82*  +  •••  +  “  0  (4.3-14) 

nuat  be  satisfied.  With  reference  to  (4.1-11),  (4.3-14)  can  be 

expanded  Into  the  following  sec  of  equations: 

j^l  jdO, 

8^+82*  +...+e  -0 

Sj+S^®  +...+e  "0 


8 


1 


These  equations  are  equivalent  to  Che  single  polynomial  equation 

8^  +  822  +  ...  +  2**  -  0  (4.3-15) 

where  the  zeros  of  (4.3-15)  are  given  by 

•  e  ;  k  -  1,2,...,  d. 

Thus,  Prony's  algebraic  equation  (4.3-8)  can  also  be  derived  by  studying 
Che  dependence  of  the  sec  of  vectors  j^,  1«1,2, . . .  ,d-fl .  As  shown  the 
DOA's  can  be  obtained  from  Its  zeros. 

We  conclude  that  Prony's  algorithm  and  the  moving  window  algorithm 
are  related  In  Che  sense  that  each  utilizes  Independence  of  a  set  of 
vectors  formed  from  the  vectors  ;  1*1,2, ... ,d't‘l. 
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4.4  Pisarenko 'a  Algorithm 


In  this  section  another  method.  Pisarenko's  algorithm  [47],  is 
presented  in  the  context  of  array  signal  processing.  Pisarenko's 
algorithm  will  also  be  shown  to  be  related  to  the  moving  rectangular 
window. 

In  Prony's  method  one  is  Interested  In  obtaining  the  roots  of 
polynomial  (4.3-8).  Since  we  are  Interested  In  the  roots  of  this 
polynomial,  the  direction  of  the  vector  ^  where 

f-ISl-Bj.  ... 

la  the  Important  criteria  and  not  Its  magnitude.  In  the  Prony  method 
we  had  set  to  I  which  constrains  the  ^  vector  to  be  on  the  hyper¬ 
plane  ■  1. 

The  goal  of  Pisarenko's  method  is  to  constrain  the  tip  of  the  ^  vector 
to  be  on  a  hypersphere.  Then  no  constraints  are  placed  in  any  direction. 

A  unit  radius  is  chosen  for  the  hypersphere  such  that 

d+l 

I  -  1.  (4.4-1) 

1-1 

Without  constraining  Co  1.  Pisarenko  solves  the  set  of  equations 

(4.3-9).  Given 


(4.4-4) 


(4.4-2)  Is  nultlplsd  by  to  obtain 

E®EB  -  0. 


Lat 


(4.4-5) 


Using  (4.4-5),  (4.4-4)  baconas 


VB  -  0.  (4.4-6) 

Noct  that  V  Is  a  (d'«-l)x(d-fl)  matrix  %rtiosa  columns  span  a  d-dlmansional 
space.  Tha  smallest  eigenvalue  of  V  must  therefore  be  equal  to  zero. 
Solving  (4.4-6)  for  £  is  equivalent  to  solving  for  tha  eiganvactor 
corresponding  to  the  smallest  (zero)  eigenvalue.  In  particular. 


VB 


^smallest  — 


•  0. 


(4.4-7) 


Based  on  this  observation,  Pisarenko  proposed  to  obtain  B  by  comput¬ 
ing  the  eigenvector  of  V  corresponding  to  the  smallest  eigenvalue. 

Once  this  eigenvector,  which  satisfies  (4.4-1),  is  obtained,  polynomial 
(4.3-8)  Is  formed  and  its  zeros  are  computed.  These  zeros, 

e  ;  i«l.,2 . d,  contain  tha  necessary  Information  to  astimate  the 

angular  position  of  tha  d  sources. 

However,  in  practice,  E^  is  not  equal  to  zero  but  to  some  residual 
£.  H.  J.  Price  [47J  solved 

E  B  •  £  (4.4-8) 

for  £  by  minimizing  the  sums  of  squares  of  the  residuals,  subject  to 
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constraint  (4,4-1).  He  showed  that  the  vector  which  ainlaizes  the 
squares  of  the  residual,  is  the  eigenvector  of  V  corresponding  to  the 
eigenvalue. 

The  relationship  between  Pisarenko's  method  and  the  moving  rec¬ 
tangular  window  is  demonstrated.  Let  be  the  column  of  V.  (4.4-6) 
can  then  be  rewritten  as 


Since  the  vectors  V^;  J«L»2.....  d+1.  are  dependent,  then  there  oust 
exist  a  set  of  constants  ft  0  such  that  (4.4-9)  is  verified.  Thus 
Pisarenko's  algorithm  examines  the  dependence  of  the  set 

‘ii-  h . 


(4.4-10) 


where 


Zj  ■  . 


j  ■  1.2,. . .  ,d4>l. 


(4.4-11) 


Recall  chat  ;  J  <■  1.2.....  d'I'l,  are  the  vectors  used  In  the  moving 
rectangular  window  asthod.  Polynomial  (4,3-S)  can  also  be  obtained 
by  examining  the  dependence  of  the  sec  (4,4-10).  Making  use  of  (4,1-13). 
Che  ij^  entry  of  V  or  can  be  expressed  as 


-  S> 

.  s“  e’  e  s 


Uslog  (4.4-12) >  can  be  written  as 


e“e 

^HCcDeH^  ^(1-1) 


Define  the  matrix  R  to  be 


* 

where  (*)  denotes  the  complex  conjugate.  The  vector  ^  can  be  expressed 
in  terms  of  R  as 
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-  a  E®  E  s.  (4.4-13) 

Utilizing  (4.4-13)  in  (4.4-9) t  w«  obtain 

3^  a  E®  E  S  +  $2  R  E®  E  «  S  +  . . . 

+  6 .  a  E®  E  s  +  3.^,  a  E®  E  -  0  (4.4-14) 

a  —  d+1  ** 

or 

a  E®E  le^I  +  82*  +...+  1  -  0.  (4.4-15) 

For  this  equality  to  hold  for  an  arbitrary  choice  of  £  It  must  be  true 
Chat 

3^1+62*+...+  “  0-  (4.4-16) 

Using  (4.1-11),  this  matrix  equation  can  be  expanded  into 

jiJ,  jdij), 

3^+82.  U...+3^^^e  ^-0 

S,  +  6.,  e  ^  +  ...  ■*•  6,_^,  e  ^  -  0  (4.4-17) 

i  z  d^i 


8.  i-  8. 


+  8 


jd(t>. 


d+1 


■  0. 


This  sec  of  equation  Is  equivalent  to  the  single  polynomial  equation. 

8^  +  822+  ...  +  8^^^  •  0  (4.4-18) 

vhere  the  zeros  of  (4.4-18)  are  given  by 
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r 


j4>. 

^  ;  i-1.2 . d.  (4.4-19) 

Indeed  polynomial  (4.4-18)  which  wee  derived  by  analyzing  che  dependence 
of  Che  aec  (4.4-10)  is  che  same  as  che  polynomial  given  in  (4.3-8). 

One  can  chen  conclude  chac  Pisarenko's  algorichm  and  che  moving 
window  algorichm  are  relaced  in  che  sense  chac  each  ucllizes  che  depen¬ 
dence/Independence  of  a  aec  of  veccors  formed  from  che  same  veccors 

X.;  i-1,2 . d+1. 

“1 
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4.S  CoapMratlv  PTfor— nee 

la  this  ssetion  chs  compsrstlvs  perfornsnee  of  ths  moving 
window,  Prony,  and  Pisarenko  methods  is  evaluated  by  means  of 
computer  simulation. 

The  model  used  for  the  slmulaClon  consists  of  two  plans  waves 
(d«2)  arising  from  two  coherent  point  sources  which  are  incident  on 
a  linear  array  consisting  of  eight  (ap8)  equally  spaced  sensor 
elements.  The  sources  are  assumed  to  be  located  at  8^  ■  18*  and 
0^  *  22*.  The  angular  separation,  A6  ■  4*,  la  lass  than  one  fourth 
of  the  array's  Rayleigh  angular  resolution  which,  for  the  given  array. 
Is  about  -  0.28  radians  >  16.3  degrees.  Assume  M  snapshots 

are  available  for  processing  where  the  snapshot  la  given  by 

^  •  {y ^ » • • • » 

For  a  given  snapshot,  the  measurement  at  the  1~  sensor  is  composed 
of  signal  and  noise  components  as  defined  In  (4.1-3).  Hence, 

y.  (t,0)  -  I  8.(t)  exp  D(l-l)  sin  9.  }  +  n  (t) 

^  k-l 

-  J  \  (1-0  slo  +  n^(t) 

k*l 

where  a^(t)  is  generated  in  the  computer  simulation  as  a  zero  mean, 
unit  variance,  white  complex  Gaussian  noise.  The  signal  portion  is 
generated  using  the  expression 
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d  n 

^  j  (i-1)  •!» 

k^l 


wfa«r«  0  <■  X/2,  d  ■  2,  dj^  •  18*,  62  "  m  ■  8.  It  is  atsuaed 

chat  cha  N  anapshots  ara  cakan  at  a  rata  aufflciantly  faat  auch  chat 
a|^(c);  k*l,2  ramaln  approxtaacaly  cooataac  ovar  cha  M  anapahoca. 
Aaaualng  Cha  2  aourcaa  to  ba  of  aqual  powar  with  aj^(c)  ■  a2(t)  ■  a, 
cha  algaal-co'^laa  ratio  SMR  la  daflnad  to  ba 

P 

SNR  - 


whara  la  cha  aignal  powar  at  aanaor  1  and  P^  la  cha  oolaa  powar 


sm 


2|a|^ 


Tha  caaaa  cenaldarad  la  cha  almulaclon  ara  cabulacad  in  Table  1. 


SNE(dB) 

Isi 

30 

22.36 

23 

12.57 

20 

7.07 

10 

2.24 

TABLE  1 


Ualng  cha  aourca  and  racalvar  aodal  aa  daacrlbad  above,  Cha  moving 
window,  Prony,  and  Plaaranko  machoda  vara  uaad  to  aaciaace  6^  and  62 
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Quantlcatlv*  results  of  the  three  aethods  ere  given  In  Tables  2 
end  3.  10  Monte  Carlo  runs  of  500  snapshots  each  were  performed  for 

each  SMR  In  the  range  of  10  -  30  dB.  It  can  be  observed  from  Table  2 
that,  at  large  enough  SHR.  the  angular  positions  of  the  two  sources 
are  nicely  resolved.  Ac  lower  SNR.  t%ro  sources  are  still  resolved 
but  Che  bias  and  variance  of  Che  esclmaces  of  6^^  and  62  ere  larger. 
Observe  that  the  moving  window  using  generalized  eigenvalues  performs 
better  chan  the  Gram  approach. 

An  unexpected  result  Is  Chat  Che  Prony  method  and  the  moving 
window  method  using  the  generalized  elgenvelue  approach  have  Identical 
performance.  This  Is  explained  below  by  showing  chat  Identical  equa¬ 
tions  are  solved  In  both  techniques.  Let 


h 


*2*  ***  •  *6^ 


^2  ■  t»2*  *3*  *  *7^ 

^  -  {Xj,  X^,  ...  ,  Xg}. 


Given  d  «  2  and  the  vectors  and  the  Prony  method  leads  to 

Che  equation 


(4.5-1) 


The  least  squares  solution  of  the  ebove  set  of  llneer  equations 
Is  given  by 
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Recall  that  tha  linear  product  of  two  vectors  and  X^  is  defined  to  be 

*2^  “52^1  “  *21*  (^.5-3) 

Using  (4.5-3) •  (4.5-2)  can  be  rewritten  as 


Utilizing  CraoMr's  rule,  and  a2  are  found  to  be 

^  ,  ^13  ^22  ~  ^23  ^12 
'  *12  *21  -  *22  *11 

and 

^23  *11  ~  *13  *21  . 

®2  "  X  Y  _  X  X 

i  *21  *12  22  11 

This  results  in  the  quadratic  equation 

Oj  \^  +  02^  +  I  -  0. 

Substituting  for  and  (X2  end  sioplifying  yields  the  equation 


(4.5-4) 


^*13*22  "  *23*12^^^  ^*23*11  '  *l3*2p^  ”*'  ^^21*12  "  *22*11^  “  ^^-5-5) 

lh€  •quaclon  In  thn  novlng  window  approach  using  ganarallzad 
algan'valuas  la  now  davalopad.  For  this,  the  matrix  pencil 


t” 

t  t 

M  -  XN  - 

1  i 

X,  X- 

-X 

1  1 

X-  X, 

—1  —2 

1  1 

—2  —3 

1  1 

I  i 

▼  Jr 

la  first  fornad.  M-XM  la  then  preoultlplled  by  yielding 


M®M-XM®N  • 

*11 

*12“ 

-X 

“*12 

*13“ 

*21 

*22 

_*22 

*23 

The  generalized  eigenvalues  are  the  solution  to 


1m®h-xm®n1  - 


*11  ‘  ^*12 


*21  '  ^22 


*12  •  '*13 


*22  "  '*23 


(*12*22  -  *23*12^^  ^*23*11  ”  *13*21^^  ^  ^*21*12  ”  *22*11^  " 


Mote  that  (4.5-5)  and  (4.5-6)  are  Identical  equations.  Consequently, 
the  two  algorithms  produce  the  same  results. 

From  Table  2  It  Is  seen  that  the  bias  and  variance  of  all  three 
methods  is  large  at  10  dB.  To  investigate  whether  9^  and  can  be 
better  estimated  using  additional  snapshots,  10  Monte  Carlo  runs  were 
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p«rfota«d  with  th*  nuab«r  of  snapshots  Incrsassd  from  SOO  to  SOOO 
whlla  SHE  was  aalntalnad  at  10  dB.  Tha  rasults  are  shorn  la  Table  4. 

As  tha  numbar  of  snapshots  Is  Incraasad,  tha  accuracy  of  cha  astimates 
for  both  anglas  also  Ineraasas.  Corraspondlngly ,  cha  varianca  given 
In  Table  S  dacraasas  as  cha  numbar  of  snapshots  incraasa.  Tha  reason 
for  this  Improvamanc  is  the  wall  lcno%m  fact  chat  cha  variance  of  the 
unbiased  aatlaacor 

i  (t^.9) 

dacraasas  as  1/M  thus  yielding  batter  results  as  N  Increases. 

Tables  6  and  7  show  cha  parformanca  of  chase  methods  as  A6  is 
successively  dacraasad  to  3*.  2*  and  1*.  For  sufficiently  large  SMR, 
tha  resolution  capabilities  of  these  methods  compare  wall.  Moca 
chat  whan  A9  Is  !*•  all  methods  fall  to  resolve  cha  2  sources.  In 
this  case,  only  one  source  Is  observed,  and  tha  astlmatad  location  of 
this  source  has  a  large  bias  and  varianca. 

Tables  8  and  9  give  Cha  results  of  the  three  methods  as  cha  sources 
are  kept  4”  apart  and  are  suvad  from  broadside  to  endflre  of  the  array. 
Toward  cha  endflre  of  tha  array,  DOA  estimation  becomes  worse.  This 
Is  because  cha  gain  of  cha  array  Is  greatly  reduced  at  endflre. 

Computer  simulations  ware  also  carried  out  to  Investigate  the 
performance  of  the  moving  window  as  given  In  section  4.1-2  where  the 
signals  are  assiimed  to  be  random  with  zero  mean.  Its  performance  Is 
compared  to  chat  of  the  ESPRIT  algorithm  under  similar  conditions.  The 
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IQQ 


Table  4:  Mean  of  Oj  and  0^  as  a  function 

of  tbe  nuiwi>er  of  anapsliots  (10  Monte 
Carlo  runs,  10  dB). 
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Table  7:  variance  of  (Ij  and  0^  ^  function 

of  AO  (500  Hnapshots,  10  Monte  Carlo 
runs,  30  dD) . 


PIOMT 


PISA&EMKO 


MOVING  WINDOW 


_ _ 

GRAM 

GEN. 

EIGEN. 

®1‘®2 

5l 

*2 

5. 

*2 

^2 

^2 

B 

.02 

4.05 

.05 

4.02 

-.30 

4.39 

mm 

8.12 

7.89 

11.90 

7.93 

11.87 

7.56 

12.24 

7.89 

11.90 

18.22 

18.02 

22.01 

18.06 

21.97 

17.62 

22.43 

18.02 

22.01 

28.32 

27.71 

31.80 

m 

27.07 

32.48 

27.71 

31.80 

36.40 

35.97 

39.99 

36.02 

39.94 

34.86 

41.18 

35.97 

39.99 

40.44 

39.97 

43.91 

40.05 

43.82 

38.50 

45.52 

39.97 

43.91 

46.50 

45.39 

49.71 

45.41 

■1 

mmmm 

43.12 

52.32 

45.39 

49.71 

54.58 

53.79 

56.47 

m 

48.45 

63.35 

53.79 

56.47 

62,66 

61.53 

65.51 

61.55 

65. 6( 

-10.99 

64.33 

61.53 

65.51 

70.74 

7.06 

74.18 

24.39 

74.0: 

-52.16 

54.61 

7.06 

74.18 

74.78 

75.85 

-8.32 

75.8: 

-45.97 

52.11 

-4.68 

75.85 

82,86 

-7.65 

83.67 

■M 

■B 

■IP 

-7.65 

83.67 

86,90 

3.68 

87.17 

mm 

liBi 

-35.93 

pp 

■Kaii 

3.68 

m 

Tabl«  8:  Maan  of  6^^  and  82  aa  cha  aourcaa  nova  from  broadside  to 
andfira  (500  snapshots.  10  Monca  Carlo  runs,  30  dB) . 
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■ 

PISARENKO 

MOVING 

f 

WINDOW  1 

GRAM 

GEN. 

EIGEN. 

B 

IH 

m 

.08 

.08 

.08 

.11 

.10 

.08 

.08 

8,12 

.03 

.03 

.04 

.03 

.06 

.04 

.03 

.03 

18,22 

.10 

.09 

.11 

.09 

.09 

.10 

.10 

.09 

28,32 

.07 

.05 

.07 

.06 

.16 

.09 

.07 

.05 

36,40 

.13 

.14 

.16 

.14 

.34 

.53 

.13 

.14 

40,44 

.11 

.28 

.13 

.36 

.47 

1.02 

.11 

.28 

46,50 

1.03 

.47 

.90 

.55 

1.75 

2.47 

1.03 

.47 

54,58 

3.01 

.29 

1.80 

.26 

5.18 

9.29 

3.01 

.29 

62,66 

6.07 

4.21 

4.44 

3.75 

3852.61 

237.12 

6.07 

4.21 

70,74 

3025.31 

26.68 

3591.77 

11.28 

505.57 

95.99 

3025.31 

26.68 

74,78 

1277.36 

.0006 

2005.98 

.011 

821.72 

504.60 

1277.36 

.0006 

82,86 

318.77 

.0006 

835.92 

.001 

1454.45 

1508.38 

318.77 

.0006 

86,90 

447.56 

.003 

1166.39 

.005 

1444.15 

1557.80 

447.56 

.003 

Table  9:  variance  of  6^  and  @2  **  sources  oove  froa  broadside  to 
endfire  (SOO  snapshots,  10  Honce  Carlo  runs,  30  dB) . 
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■odd  used  for  sloulstion  is  the  ssas  ss  the  ons  described  above 
except  Chet  Che  sources  ere  now  asde  incoherent.  The  reason  for 
this  is  Che  ESPRIT  algorithu  fails  for  coherent  sources. 

Making  Che  sources  incoherent  gives  us  a  fair  coaperison  of  the 
■ovlng  window  and  ESPRIT.  The  signals  are  generated  by  using  the 
expression. 

2 

J  “P  (i“l)  sin  0  }  exp  {j  ol  }, 

k*l 

and  Che  noise  is  generated  as  before.  Two  cases  were  investigated  in  this 
sioulacion.  In  case  1,  is  kept  constant  and  and  02  arc  inde¬ 

pendent  randoa  phase  angles  uniformly  distributed  in  the  interval 
In  case  2,  s^(c)  is  randoa  with  a  Rayleigh  distribution  independent  of 
and  02 t  while  o^  and  O2  are  still  independent  randoa  phase  angles 
uniforaly  distributed  in  the  interval  [-r,  r].  Tables  10  and  11  show 
Che  results  for  case  1,  and  Tables  12  and  13  show  chose  of  case  2. 

The  aacrlces  M  and  N  for  the  aovlng  window  are  forasd  as  explained 
in  section  4.1-2.  For  the  ESPRIT  algoritha,  the  first  subarray  is 
formed  of  the  first,  third,  fifth  and  seventh  sensors  while  the  second 
subarray  is  formed  of  the  second,  fourth,  sixth  and  eighth  sensors  of 
Che  linear  array.  These  two  subarrays  are  used  to  form  the  covariance 
matrices  ESPRIT  calls  for.  For  every  Honte  Carlo  run,  Che  matrix 
entries  for  both  techniques  were  computed  from  100  snapshots.  SO 
Monte  Carlo  runs  were  performed  using  independent  data  sets. 
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On*  can  conclude  froa  table  10  -  13  chat  Che  aovlng  window 
coaparea  favorably  to  ESP&IT.  The  bias  Is  slightly  saaller  for 
Che  aovlng  window  technique  while  Its  variance  is  substantially  saaller. 
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CHAPTER  S 


CONCLUSION  AND  SUGGESTIONS 
FOR  FUTURE  RESEARCH 


5.1  Concltt«loti 

D«t*tmlnatlon  of  th«  angular  locations  of  d  sourcas  using  an 
array  of  a  sansors  was  cha  aaln  concanx  of  this  ros&arcb..  Tha 

signals  ganaratad  by  cha  d  sourcas  ara  assuaad  to  ba  narrowband. 

Until  racantly  cha  asthods  davalopad  to  daal  with  this  problaa 
wars  classlflad  as  saarch  procaduras  bacauaa  althar 

1)  cha  algorltha  solvas  a  eonstralnad  opClad.zatlon  problaa  for 
aach  dlracclon  of  look, 

2)  a  baaa  Is  foxaad  and  Its  anargy  coaputsd  for  aach  dlracclon 
of  lookt  or 

3)  a  spatial  corralatlon  aacrtz  Is  foraad.  Tha  array  aanlfold 
Is  chan  saarchad  for  cha  yaluaa  of  a (6)  which  alnlalza  a  prsdaflnad 
axprasslon. 

Although  thasa  saarch  tachnlquas  nay  have  superresoluclon  capabll- 
Idas,  chay  ara,  navarthalass ,  coapuatlonally  vary  coaq>lax.  Bacause 
of  this  conplaxlty,  nonsaarch  procaduras  hava  baan  proposad.  Thasa 
algorlthas  haws  cha  following  advancagaa  ovar  saarch  procaduras: 

1)  Thay  ara  conputatlonally  lass  conplax  bacausa  a  saarch  pro- 
cadura  Is  not  naadad. 

2)  Thay  do  not  raqulra  knowladga  of  alanant  charactarlstlcs . 

3)  Thay  do  not  raqulra  a  calibration  of  tha  array.  This 
conplacaly  allialnacas  naad  for  storaga  of  cha  array  aanlfold  which 
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can  b«  vary  larg*  for  aulcldlaanalonal  problaaa. 

4)  Thoy  nay  not  raqulra  kaofwladga  of  tha  array  gaooiacry,  as  Is 
cha  caaa  with  ESPRIT. 

A  ganarallzad  foraulaclon  la  propoaad  for  thaaa  nonaaarch  proeaduras. 
Thla  foraulaclon  conslaca  of  forming  a  natrlx  pancll  M-XN  and  conputlng 
Ics  rank  raduclog  nuabara.  lha  rank  raduelng  miabara  ara  cha  ganarallzad 
alganvaluas  for  cha  caaa  of  a  aquara  nacrlx  pancll.  For  cha  caaa  of 
a  raccangular  oacrlx  pancll t  cha  rank  raduelng  nuabara  can  ba  obcalnad 
by  alchar  using  cha  Granalan  approach  or  by  cransforalng  cha  problaa  Inco 
a  ganarallzad  alganvalua  problaa.  Tha  pancll  chaoraa  glvan  In  chapear 
2  aacabllshaa  cha  ralaclonshlp  baevaan  cha  rank  raduelng  valuas  of  X 
and  a  funeclonal  fom  f(4^)  vhich  la  a  nonllnaar  funcClon  of  Cha 
angular  poalclon  of  Cha  1~  aourea;  l*l>2>...»d.  Tha  fora  of  f(d^)  la 
dacaralnad  by  cha  oparacors  appllad  co  cha  aaaauraaancs . 

Ihraa  dlffaranc  oparacors  ara  prasancad  In  chls  dlasarcaclon:  cha 
phaaa  daisy  oparacor,  cha  suaaaclon  oparacor,  and  cha  aovlng  window 
oparacor.  All  Chraa  aachods  ara  analyzad  and  fomilacad  In  earns  of 
cha  ganarallzad  foraulaclon.  Tha  aacrlz  dacoaposlclon  raquirad  by 
cha  pancll  chaoraa  haps  co  axplaln  why  ESPRIT  falls  for  coharanc  sig¬ 
nals.  In  parclculsTt  cha  rank  raqulraasnc  falls  whan  cha  signals  are 
coharanc.  Tha  moving  window  and  cha  suaaaclon  oparacor  weta  sho«m  noc 
co  vlolaca  any  of  cha  pancll  chaoraa  raqulraaancs  and  chus  do  noc  fall 
for  coharanc  signals. 

Tha  flazlblllcy  of  cha  aacrlx  pancll  approach  allowad  us  co  forau- 
laca  cha  aovlng  window  oparacor  for  cha  casa  of  dacamlnlsclc  signals  as 
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wall  as  for  tha  zaro  oaan  random  caaa.  Two  wall  known  aathoda  In  tha 


flald  of  ayatam  Idantlflcation*  naaaly  Prony'a  and  Plaaranko's  nachoda, 
ara  appllad  to  tha  problam  of  diractlon  finding.  Thalr  ralatlonahlp 
to  tha  aovlng  window  tachnlqua  la  daannatratad .  All  thraa  tachnlquas 
aztraet  tha  Information  about  tha  angular  loeatlona  of  tha  d  algnala 
by  axamlnlng  tha  dapandanea/lndapandanca  of  a  aat  of  vactors.  Each 
tachnlqua  darlvaa  tha  aat  of  vaetora  aa  a  tranaforMtlon  of  tha  aama 
(d  -f  1)  vaetora. 

Tha  nonaaareh  aathoda  praaanead  la  thla  research  ara  capable 
of  making  high  raaolutlon  OOA  aaelaatlon.  lhaaa  aathoda  hava  great 
potential  baeauaa  of  thalr  computational  alapllclty.  They  ara  aaay 
and  cheap  to  laplamant. 

Aa  one  davalopa  an  Idea,  there  ara  alwaya  new  quaatlona  and  prob- 
lama  that  arlaa  In  tha  procaaa.  In  tha  following  aaetlon  auggaatlona 
ara  made  for  future  raaaarch. 

S.2  Suaaaatlona  for  Future  Raaaarch 

The  genarallzad  approach  la,  aa  Ita  name  Indlcataa,  a  vary  general 
concapt.  Tha  rank  reducing  numbara,  X2,...,  X^,  of  the  matrix 
pencil  H-XN  ara  related  to  tha  dlrectlona  of  arrival  thru  a  functional 
form  that  dapanda  on  the  operator  applied  to  the  meaaurementa.  The 
poaalblllty  of  ualng  operatora  other  than  thoae  dlacuaaad  in  this 
diaaartatlon  raaalna  to  be  explored. 

Tha  formulation  of  different  aathoda  In  terna  of  a  conon  frame¬ 
work  haa  the  advantage  that  their  performance  can  be  compared  within 


th«  ganarallzad  approach.  For  axaapla.  tupposa  cm  oachoda  arc  fomu- 
latad  in  eons  of  cho  aacrlx  ponell.  Procosslng  of  cha  aacrlz  pancll 
la  Idantlcal  for  both  MChods.  Bowavar*  cha  ancrlas  in  cha  pancll  vill 
dlffar.  k  acaelatical  analysis  of  chasa  ancrlas  can  chon  ba  carriad 
ouc.  Iho  naehod  whoaa  ancrlas  hava  saallar  bias  and  saallar  varlanca 
will  parfon  baetar. 

Ua  know  fron  cha  dasign  of  FIX.  fllcars  chac  dlffaranc  trlndovs  hava 
dlffaranc  characcarlsclcs .  Tha  window  usad  In  chapcar  4  is  a  racCangtilar 
window.  Ic  can  vary  aaslly  ba  shown  chac  any  shapa  window  would  mrk. 

Tha  quaaclon  Is*  "Bow  do  dlffarancly  shapad  windows  affacc  cha  par- 
forwanca  of  Cha  aovlng  window  oparacor?" 

Ic  was  shown  In  sacclon  4.1-3  chac  prof 11 caring  can  ba  Inciudad  Inco 
Chasa  nonsaarch  procaduras  wlchouc  dlscurblng  cha  nacrlx  dacoaposlCion 
raqulrad  by  cha  pancll  chaoran.  An  laporcanc  Issuo  concarnlng  pra- 
f Hearing  la  cha  dasign  of  a  sulcabla  fllcar.  If  ona  know  anough  abouc 
cha  nolsa  and  cha  daslrad  signal*  ona  muld  ba  abla  co  build  a  filcar 
CO  laprova  cha  signal -co-nolsa  raclo  and,  charoforo*  laprova  cha  par- 
foraanca  of  cha  cachnlqua. 

Tha  signals  In  dlracclon  finding*  spaccral  asclaaclon,  syscaa  Idancl- 
flcaclon*  and  adaptlva  arrays  can  all  ba  aodalad  as  a  sua  of  axponanclals 
Conaaquontly*  afforc  should  ba  davoead  co  applying  cha  rasulca  of  chis 
dlaartaclon  co  chosa  oehor  araas. 

Anochar  Issua  chac  should  ba  glvan  consldaraclon  Is  sansor  coupling. 
Tha  affacc  of  sansor  coupling  has  noc  baan  dlscussad  In  cha  opan  licar- 
acuca  on  dlracclon  finding.  Invasclgaclon  Inco  chis  problan  alghc  help 
laprova  axlsclng  algorlehaa. 


Tha  gtttMrallzad  fraatvork  d«vclop«d  In  this  research  relics 
on  the  coaputatlon  of  the  rank  reducing  nimbers  of  a  matrix  pencil 
M-XN.  An  efficient  algorithm  [48]  baaed  on  the  Schur  decomposition 
has  been  developed  for  the  case  of  a  square  matrix  pencil.  However, 
for  the  case  of  a  rectangular  pencil,  no  efficient  algorithm  has  been 
shown  for  explicitly  computing  the  coefficients  of  the  polynomial 
P(X)  •  det  (G)  *0.  Preliminary  work  suggests  that  such  an  approach 
la  possible  [42.  47]. 
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